Preparatory steps

Load packages

required_packages = c("confintr", "tidyr","ggplot2", "effectsize", "easystats", "Jmisc", "sjmisc", "plyr", "lme4", "lmerTest", "emmeans",  "dplyr", "ggpubr", "purrr", "broom",  "plotrix", "VGAM", "reshape2", "PupillometryR")
#detach(package:lmerTest, unload=TRUE)
#library(scales)
#library(ggplot2)
lapply(required_packages, library, character.only = TRUE)
## 
## Attaching package: 'effectsize'
## The following object is masked from 'package:confintr':
## 
##     oddsratio
## # Attaching packages: easystats 0.6.0
## ✔ bayestestR  0.13.0   ✔ correlation 0.8.3 
## ✔ datawizard  0.6.5    ✔ insight     0.19.0
## ✔ modelbased  0.8.6    ✔ performance 0.10.2
## ✔ parameters  0.20.2   ✔ report      0.5.6 
## ✔ see         0.7.4
## 
## Attaching package: 'Jmisc'
## The following object is masked from 'package:parameters':
## 
##     demean
## The following object is masked from 'package:datawizard':
## 
##     demean
## The following object is masked from 'package:ggplot2':
## 
##     %+%
## 
## Attaching package: 'sjmisc'
## The following objects are masked from 'package:datawizard':
## 
##     center, empty_rows, remove_empty_rows, reshape_longer, to_factor,
##     to_numeric
## The following object is masked from 'package:tidyr':
## 
##     replace_na
## 
## Attaching package: 'plyr'
## The following object is masked from 'package:here':
## 
##     here
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:Jmisc':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
## 
##     mutate
## 
## Attaching package: 'purrr'
## The following object is masked from 'package:plyr':
## 
##     compact
## The following object is masked from 'package:sjmisc':
## 
##     is_empty
## 
## Attaching package: 'plotrix'
## The following object is masked from 'package:datawizard':
## 
##     rescale
## Loading required package: stats4
## Loading required package: splines
## 
## Attaching package: 'VGAM'
## The following object is masked from 'package:tidyr':
## 
##     fill
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
## Loading required package: rlang
## 
## Attaching package: 'rlang'
## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice
## The following object is masked from 'package:sjmisc':
## 
##     is_empty
## [[1]]
##  [1] "confintr"  "knitr"     "here"      "pacman"    "stats"     "graphics" 
##  [7] "grDevices" "datasets"  "utils"     "methods"   "base"     
## 
## [[2]]
##  [1] "tidyr"     "confintr"  "knitr"     "here"      "pacman"    "stats"    
##  [7] "graphics"  "grDevices" "datasets"  "utils"     "methods"   "base"     
## 
## [[3]]
##  [1] "ggplot2"   "tidyr"     "confintr"  "knitr"     "here"      "pacman"   
##  [7] "stats"     "graphics"  "grDevices" "datasets"  "utils"     "methods"  
## [13] "base"     
## 
## [[4]]
##  [1] "effectsize" "ggplot2"    "tidyr"      "confintr"   "knitr"     
##  [6] "here"       "pacman"     "stats"      "graphics"   "grDevices" 
## [11] "datasets"   "utils"      "methods"    "base"      
## 
## [[5]]
##  [1] "see"         "report"      "parameters"  "performance" "modelbased" 
##  [6] "insight"     "datawizard"  "correlation" "bayestestR"  "easystats"  
## [11] "effectsize"  "ggplot2"     "tidyr"       "confintr"    "knitr"      
## [16] "here"        "pacman"      "stats"       "graphics"    "grDevices"  
## [21] "datasets"    "utils"       "methods"     "base"       
## 
## [[6]]
##  [1] "Jmisc"       "see"         "report"      "parameters"  "performance"
##  [6] "modelbased"  "insight"     "datawizard"  "correlation" "bayestestR" 
## [11] "easystats"   "effectsize"  "ggplot2"     "tidyr"       "confintr"   
## [16] "knitr"       "here"        "pacman"      "stats"       "graphics"   
## [21] "grDevices"   "datasets"    "utils"       "methods"     "base"       
## 
## [[7]]
##  [1] "sjmisc"      "Jmisc"       "see"         "report"      "parameters" 
##  [6] "performance" "modelbased"  "insight"     "datawizard"  "correlation"
## [11] "bayestestR"  "easystats"   "effectsize"  "ggplot2"     "tidyr"      
## [16] "confintr"    "knitr"       "here"        "pacman"      "stats"      
## [21] "graphics"    "grDevices"   "datasets"    "utils"       "methods"    
## [26] "base"       
## 
## [[8]]
##  [1] "plyr"        "sjmisc"      "Jmisc"       "see"         "report"     
##  [6] "parameters"  "performance" "modelbased"  "insight"     "datawizard" 
## [11] "correlation" "bayestestR"  "easystats"   "effectsize"  "ggplot2"    
## [16] "tidyr"       "confintr"    "knitr"       "here"        "pacman"     
## [21] "stats"       "graphics"    "grDevices"   "datasets"    "utils"      
## [26] "methods"     "base"       
## 
## [[9]]
##  [1] "lme4"        "Matrix"      "plyr"        "sjmisc"      "Jmisc"      
##  [6] "see"         "report"      "parameters"  "performance" "modelbased" 
## [11] "insight"     "datawizard"  "correlation" "bayestestR"  "easystats"  
## [16] "effectsize"  "ggplot2"     "tidyr"       "confintr"    "knitr"      
## [21] "here"        "pacman"      "stats"       "graphics"    "grDevices"  
## [26] "datasets"    "utils"       "methods"     "base"       
## 
## [[10]]
##  [1] "lmerTest"    "lme4"        "Matrix"      "plyr"        "sjmisc"     
##  [6] "Jmisc"       "see"         "report"      "parameters"  "performance"
## [11] "modelbased"  "insight"     "datawizard"  "correlation" "bayestestR" 
## [16] "easystats"   "effectsize"  "ggplot2"     "tidyr"       "confintr"   
## [21] "knitr"       "here"        "pacman"      "stats"       "graphics"   
## [26] "grDevices"   "datasets"    "utils"       "methods"     "base"       
## 
## [[11]]
##  [1] "emmeans"     "lmerTest"    "lme4"        "Matrix"      "plyr"       
##  [6] "sjmisc"      "Jmisc"       "see"         "report"      "parameters" 
## [11] "performance" "modelbased"  "insight"     "datawizard"  "correlation"
## [16] "bayestestR"  "easystats"   "effectsize"  "ggplot2"     "tidyr"      
## [21] "confintr"    "knitr"       "here"        "pacman"      "stats"      
## [26] "graphics"    "grDevices"   "datasets"    "utils"       "methods"    
## [31] "base"       
## 
## [[12]]
##  [1] "dplyr"       "emmeans"     "lmerTest"    "lme4"        "Matrix"     
##  [6] "plyr"        "sjmisc"      "Jmisc"       "see"         "report"     
## [11] "parameters"  "performance" "modelbased"  "insight"     "datawizard" 
## [16] "correlation" "bayestestR"  "easystats"   "effectsize"  "ggplot2"    
## [21] "tidyr"       "confintr"    "knitr"       "here"        "pacman"     
## [26] "stats"       "graphics"    "grDevices"   "datasets"    "utils"      
## [31] "methods"     "base"       
## 
## [[13]]
##  [1] "ggpubr"      "dplyr"       "emmeans"     "lmerTest"    "lme4"       
##  [6] "Matrix"      "plyr"        "sjmisc"      "Jmisc"       "see"        
## [11] "report"      "parameters"  "performance" "modelbased"  "insight"    
## [16] "datawizard"  "correlation" "bayestestR"  "easystats"   "effectsize" 
## [21] "ggplot2"     "tidyr"       "confintr"    "knitr"       "here"       
## [26] "pacman"      "stats"       "graphics"    "grDevices"   "datasets"   
## [31] "utils"       "methods"     "base"       
## 
## [[14]]
##  [1] "purrr"       "ggpubr"      "dplyr"       "emmeans"     "lmerTest"   
##  [6] "lme4"        "Matrix"      "plyr"        "sjmisc"      "Jmisc"      
## [11] "see"         "report"      "parameters"  "performance" "modelbased" 
## [16] "insight"     "datawizard"  "correlation" "bayestestR"  "easystats"  
## [21] "effectsize"  "ggplot2"     "tidyr"       "confintr"    "knitr"      
## [26] "here"        "pacman"      "stats"       "graphics"    "grDevices"  
## [31] "datasets"    "utils"       "methods"     "base"       
## 
## [[15]]
##  [1] "broom"       "purrr"       "ggpubr"      "dplyr"       "emmeans"    
##  [6] "lmerTest"    "lme4"        "Matrix"      "plyr"        "sjmisc"     
## [11] "Jmisc"       "see"         "report"      "parameters"  "performance"
## [16] "modelbased"  "insight"     "datawizard"  "correlation" "bayestestR" 
## [21] "easystats"   "effectsize"  "ggplot2"     "tidyr"       "confintr"   
## [26] "knitr"       "here"        "pacman"      "stats"       "graphics"   
## [31] "grDevices"   "datasets"    "utils"       "methods"     "base"       
## 
## [[16]]
##  [1] "plotrix"     "broom"       "purrr"       "ggpubr"      "dplyr"      
##  [6] "emmeans"     "lmerTest"    "lme4"        "Matrix"      "plyr"       
## [11] "sjmisc"      "Jmisc"       "see"         "report"      "parameters" 
## [16] "performance" "modelbased"  "insight"     "datawizard"  "correlation"
## [21] "bayestestR"  "easystats"   "effectsize"  "ggplot2"     "tidyr"      
## [26] "confintr"    "knitr"       "here"        "pacman"      "stats"      
## [31] "graphics"    "grDevices"   "datasets"    "utils"       "methods"    
## [36] "base"       
## 
## [[17]]
##  [1] "VGAM"        "splines"     "stats4"      "plotrix"     "broom"      
##  [6] "purrr"       "ggpubr"      "dplyr"       "emmeans"     "lmerTest"   
## [11] "lme4"        "Matrix"      "plyr"        "sjmisc"      "Jmisc"      
## [16] "see"         "report"      "parameters"  "performance" "modelbased" 
## [21] "insight"     "datawizard"  "correlation" "bayestestR"  "easystats"  
## [26] "effectsize"  "ggplot2"     "tidyr"       "confintr"    "knitr"      
## [31] "here"        "pacman"      "stats"       "graphics"    "grDevices"  
## [36] "datasets"    "utils"       "methods"     "base"       
## 
## [[18]]
##  [1] "reshape2"    "VGAM"        "splines"     "stats4"      "plotrix"    
##  [6] "broom"       "purrr"       "ggpubr"      "dplyr"       "emmeans"    
## [11] "lmerTest"    "lme4"        "Matrix"      "plyr"        "sjmisc"     
## [16] "Jmisc"       "see"         "report"      "parameters"  "performance"
## [21] "modelbased"  "insight"     "datawizard"  "correlation" "bayestestR" 
## [26] "easystats"   "effectsize"  "ggplot2"     "tidyr"       "confintr"   
## [31] "knitr"       "here"        "pacman"      "stats"       "graphics"   
## [36] "grDevices"   "datasets"    "utils"       "methods"     "base"       
## 
## [[19]]
##  [1] "PupillometryR" "rlang"         "reshape2"      "VGAM"         
##  [5] "splines"       "stats4"        "plotrix"       "broom"        
##  [9] "purrr"         "ggpubr"        "dplyr"         "emmeans"      
## [13] "lmerTest"      "lme4"          "Matrix"        "plyr"         
## [17] "sjmisc"        "Jmisc"         "see"           "report"       
## [21] "parameters"    "performance"   "modelbased"    "insight"      
## [25] "datawizard"    "correlation"   "bayestestR"    "easystats"    
## [29] "effectsize"    "ggplot2"       "tidyr"         "confintr"     
## [33] "knitr"         "here"          "pacman"        "stats"        
## [37] "graphics"      "grDevices"     "datasets"      "utils"        
## [41] "methods"       "base"
#pacman::p_load(char = required_packages)
source(here::here("utils/r_functions.R"))

Load additional tools

# load color palettes
pal <- get_colors("ond2")
pal2 <- get_colors("ond")
pal3 <- get_colors("ukr")

plot(NULL, xlim=c(0,length(pal2)), ylim=c(0,1), 
    xlab="", ylab="", xaxt="n", yaxt="n")
rect(0:(length(pal)-1), 0.7, 1:length(pal), 0.9, col=pal)
rect(0:(length(pal2)-1), 0.4, 1:length(pal2), 0.6, col=pal2)
rect(0:(length(pal3)-1), 0.1, 1:length(pal3), 0.3, col=pal3)

1. Analysis of ratings - stable cues

Stats

stable_data<-read.csv(here::here("data/stable_cues_data.csv"))
stable_data = assign_var_types(stable_data, c("trtype_str", "study_str", "id", "contingency", "ta"))

# estimate models
m1a<-lmer(prob ~ ta*contingency*trtype_str + (1|id) + (1|study_str), stable_data)
## boundary (singular) fit: see help('isSingular')
m1b<-lmer(prob ~ ta*contingency*trtype_str + (1|id) + (1+ta+contingency+trtype_str|study_str), stable_data)
## boundary (singular) fit: see help('isSingular')
# model comparison with and without random slope
anova(m1a, m1b)
## refitting model(s) with ML (instead of REML)
## Data: stable_data
## Models:
## m1a: prob ~ ta * contingency * trtype_str + (1 | id) + (1 | study_str)
## m1b: prob ~ ta * contingency * trtype_str + (1 | id) + (1 + ta + contingency + trtype_str | study_str)
##     npar     AIC      BIC logLik deviance Chisq Df Pr(>Chisq)
## m1a   15 -232.54 -175.924 131.27  -262.54                    
## m1b   29 -200.53  -91.073 129.27  -258.54     0 14          1
# winning model overview
anova(m1a)
## Type III Analysis of Variance Table with Satterthwaite's method
##                            Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
## ta                         0.0000  0.0000     1 309.818   0.0000  0.996861    
## contingency                0.0078  0.0039     2  17.227   0.1455  0.865610    
## trtype_str                11.7153 11.7153     1 308.261 435.8064 < 2.2e-16 ***
## ta:contingency             0.0270  0.0135     2 309.263   0.5023  0.605632    
## ta:trtype_str              0.1857  0.1857     1 308.261   6.9064  0.009019 ** 
## contingency:trtype_str     1.8189  0.9094     2 308.261  33.8311 5.211e-14 ***
## ta:contingency:trtype_str  0.0781  0.0391     2 308.261   1.4528  0.235503    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effectsize
effectsize::effectsize(anova(m1a), type="eta",  alternative = "two.sided")
## # Effect Size for ANOVA (Type III)
## 
## Parameter                 | Eta2 (partial) |       95% CI
## ---------------------------------------------------------
## ta                        |       5.00e-08 | [0.00, 0.00]
## contingency               |           0.02 | [0.00, 0.17]
## trtype_str                |           0.59 | [0.52, 0.64]
## ta:contingency            |       3.24e-03 | [0.00, 0.02]
## ta:trtype_str             |           0.02 | [0.00, 0.06]
## contingency:trtype_str    |           0.18 | [0.11, 0.25]
## ta:contingency:trtype_str |       9.34e-03 | [0.00, 0.04]
# statistical assumptions
## homogeneity of variance 
stable_data$residuals <- residuals(m1a)
stable_data$residuals_sq <- abs(residuals(m1a))^2
hist(stable_data$residuals_sq)

lev_m <- lm(residuals_sq ~ id, data=stable_data) #ANOVA of the squared residuals
anova(lev_m)
## Analysis of Variance Table
## 
## Response: residuals_sq
##            Df  Sum Sq   Mean Sq F value  Pr(>F)  
## id         88 0.28219 0.0032067   1.267 0.08311 .
## Residuals 233 0.58972 0.0025310                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
levRes = car::leveneTest(residuals_sq ~ id, data=stable_data)
levRes
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group  88  1.2283 0.1142
##       233
## normality of residuals
ggplot(data=stable_data, aes(x=residuals)) + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

car::qqPlot(stable_data$residuals)

## [1] 152 188
# marginals and post-hocs 
em1a1 = emmeans(m1a, specs = pairwise ~ trtype_str)
## NOTE: Results may be misleading due to involvement in interactions
em1a1$emmeans
##  trtype_str emmean    SE   df lower.CL upper.CL
##  hi-prob     0.718 0.023 1.23    0.528    0.909
##  low-prob    0.302 0.023 1.23    0.112    0.493
## 
## Results are averaged over the levels of: contingency 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
em1a1$contrasts
##  contrast               estimate     SE  df t.ratio p.value
##  (hi-prob) - (low-prob)    0.416 0.0199 242  20.876  <.0001
## 
## Results are averaged over the levels of: contingency 
## Degrees-of-freedom method: kenward-roger
em1a2 = emmeans(m1a, specs = pairwise ~ trtype_str*contingency)
## NOTE: Results may be misleading due to involvement in interactions
em1a2$emmeans
##  trtype_str contingency emmean     SE   df lower.CL upper.CL
##  hi-prob    60/40        0.608 0.0373 6.76   0.5195    0.697
##  low-prob   60/40        0.417 0.0373 6.76   0.3285    0.506
##  hi-prob    75/25        0.724 0.0185 6.24   0.6792    0.769
##  low-prob   75/25        0.306 0.0185 6.24   0.2615    0.351
##  hi-prob    90/10        0.822 0.0370 6.39   0.7327    0.911
##  low-prob   90/10        0.184 0.0370 6.39   0.0945    0.273
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
em1a2$contrasts
##  contrast                            estimate     SE    df t.ratio p.value
##  (hi-prob 60/40) - (low-prob 60/40)    0.1910 0.0387 242.1   4.932  <.0001
##  (hi-prob 60/40) - (hi-prob 75/25)    -0.1158 0.0393  49.2  -2.947  0.0522
##  (hi-prob 60/40) - (low-prob 75/25)    0.3019 0.0393  49.2   7.686  <.0001
##  (hi-prob 60/40) - (hi-prob 90/10)    -0.2135 0.0384 243.7  -5.554  <.0001
##  (hi-prob 60/40) - (low-prob 90/10)    0.4247 0.0384 243.7  11.047  <.0001
##  (low-prob 60/40) - (hi-prob 75/25)   -0.3068 0.0393  49.2  -7.810  <.0001
##  (low-prob 60/40) - (low-prob 75/25)   0.1109 0.0393  49.2   2.823  0.0705
##  (low-prob 60/40) - (hi-prob 90/10)   -0.4046 0.0384 243.7 -10.523  <.0001
##  (low-prob 60/40) - (low-prob 90/10)   0.2337 0.0384 243.7   6.078  <.0001
##  (hi-prob 75/25) - (low-prob 75/25)    0.4177 0.0247 242.1  16.881  <.0001
##  (hi-prob 75/25) - (hi-prob 90/10)    -0.0977 0.0390  47.4  -2.507  0.1426
##  (hi-prob 75/25) - (low-prob 90/10)    0.5405 0.0390  47.4  13.863  <.0001
##  (low-prob 75/25) - (hi-prob 90/10)   -0.5155 0.0390  47.4 -13.221  <.0001
##  (low-prob 75/25) - (low-prob 90/10)   0.1228 0.0390  47.4   3.149  0.0318
##  (hi-prob 90/10) - (low-prob 90/10)    0.6382 0.0381 242.1  16.732  <.0001
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 6 estimates
em1a3 = emmeans(m1a, specs = pairwise ~ trtype_str | contingency)
## NOTE: Results may be misleading due to involvement in interactions
em1a3$emmeans
## contingency = 60/40:
##  trtype_str emmean     SE   df lower.CL upper.CL
##  hi-prob     0.608 0.0373 6.76   0.5195    0.697
##  low-prob    0.417 0.0373 6.76   0.3285    0.506
## 
## contingency = 75/25:
##  trtype_str emmean     SE   df lower.CL upper.CL
##  hi-prob     0.724 0.0185 6.24   0.6792    0.769
##  low-prob    0.306 0.0185 6.24   0.2615    0.351
## 
## contingency = 90/10:
##  trtype_str emmean     SE   df lower.CL upper.CL
##  hi-prob     0.822 0.0370 6.39   0.7327    0.911
##  low-prob    0.184 0.0370 6.39   0.0945    0.273
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
em1a3$contrasts
## contingency = 60/40:
##  contrast               estimate     SE  df t.ratio p.value
##  (hi-prob) - (low-prob)    0.191 0.0387 242   4.932  <.0001
## 
## contingency = 75/25:
##  contrast               estimate     SE  df t.ratio p.value
##  (hi-prob) - (low-prob)    0.418 0.0247 242  16.881  <.0001
## 
## contingency = 90/10:
##  contrast               estimate     SE  df t.ratio p.value
##  (hi-prob) - (low-prob)    0.638 0.0381 242  16.732  <.0001
## 
## Degrees-of-freedom method: kenward-roger
a = summary(em1a3$contrasts)
effectsize::t_to_eta2(t=a$t.ratio, df_error = a$df, alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.09           | [0.03, 0.17]
## 0.54           | [0.46, 0.61]
## 0.54           | [0.46, 0.60]
# trend with anxiety
emtr<-emtrends(m1a, pairwise ~ trtype_str, var = "ta")
## NOTE: Results may be misleading due to involvement in interactions
sumtr<-summary(emtr$contrasts)
effectsize::t_to_eta2(t=sumtr$t.ratio, df_error = sumtr$df, type="eta",  alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.03           | [0.00, 0.08]
# one-way t-tests against true reinforcement (estimated by err)

df2 <- stable_data %>% 
    group_by(contingency, tabin, trtype_str) %>%
    nest() %>% 
    dplyr::mutate(tt=map(data,~t.test(.x$err))) %>%
    dplyr::mutate(tidied = map(tt, tidy)) %>% 
    unnest(tidied, .drop = T)

df2$p.value.fdr = p.adjust(df2$p.value, method = "fdr")

df2
## # A tibble: 12 × 14
## # Groups:   contingency, tabin, trtype_str [12]
##    tabin contingency trtype_str data  tt    estimate statistic p.value parameter
##    <fct> <fct>       <fct>      <lis> <lis>    <dbl>     <dbl>   <dbl>     <dbl>
##  1 lowA… 75/25       hi-prob    <tib… <hte… -0.0535    -2.62   1.18e-2        47
##  2 lowA… 75/25       low-prob   <tib… <hte…  0.0810     3.58   8.18e-4        47
##  3 hiAnx 75/25       hi-prob    <tib… <hte… -0.0344    -1.32   1.95e-1        39
##  4 hiAnx 75/25       low-prob   <tib… <hte…  0.0428     2.04   4.78e-2        39
##  5 hiAnx 60/40       hi-prob    <tib… <hte… -0.0179    -0.424  6.77e-1        17
##  6 hiAnx 60/40       low-prob   <tib… <hte…  0.0239     0.682  5.05e-1        17
##  7 hiAnx 90/10       hi-prob    <tib… <hte… -0.0723    -2.07   5.42e-2        17
##  8 hiAnx 90/10       low-prob   <tib… <hte…  0.0160     0.776  4.49e-1        17
##  9 lowA… 60/40       hi-prob    <tib… <hte… -0.00191   -0.0421 9.67e-1        17
## 10 lowA… 60/40       low-prob   <tib… <hte…  0.00247    0.0567 9.55e-1        17
## 11 lowA… 90/10       hi-prob    <tib… <hte… -0.152     -3.51   2.52e-3        18
## 12 lowA… 90/10       low-prob   <tib… <hte…  0.157      2.32   3.21e-2        18
## # … with 5 more variables: conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.value.fdr <dbl>
effectsize::t_to_eta2(t=df2$statistic, df_error = df2$parameter, alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.13           | [0.01, 0.31]
## 0.21           | [0.04, 0.40]
## 0.04           | [0.00, 0.22]
## 0.10           | [0.00, 0.30]
## 0.01           | [0.00, 0.24]
## 0.03           | [0.00, 0.29]
## 0.20           | [0.00, 0.50]
## 0.03           | [0.00, 0.31]
## 1.04e-04       | [0.00, 0.03]
## 1.89e-04       | [0.00, 0.06]
## 0.41           | [0.07, 0.65]
## 0.23           | [0.00, 0.52]
## Gather model predictions and summarize for plots below 
stable_data$model_pred <- predict(m1a) 

Plots

Stables 1: High/Low cue

stable_data_df = stable_data %>%
  group_by(id, trtype_str) %>%
  summarise_at(c("model_pred","prob"), mean, na.rm = TRUE)


g <- ggplot(data = stable_data_df, aes(y = prob, x = trtype_str, fill = trtype_str)) +
geom_hline(yintercept = 0.75, color=pal[5], linetype="longdash")  +
geom_hline(yintercept = 0.25, color=pal[7], linetype="longdash")  +
  
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2) +
geom_point(aes(y = prob, color=trtype_str), position = position_jitter(width = .15), size = 2, alpha = 0.5, show.legend = FALSE) +
geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +
  

stat_summary(geom="point", fun="mean", size=4, shape=23,  stroke=1.5, aes(y=model_pred, group=trtype_str), color="black", show.legend = FALSE) +

expand_limits(x = 3) +

scale_color_manual(values = c(pal[5], pal[7])) +
scale_fill_manual(values = c(pal[5], pal[7]), name = "Stable cue", labels = c("High-prob Cue", "Low-prob Cue")) +  
# coord_flip() +
theme_bw() +
raincloud_theme +

theme(legend.position = "right") + 
labs(y= "Mean probability", x="Session") +
  scale_x_discrete(labels=c("hi-prob" = "HC", "low-prob" =  "LC")) + 
guides(fill=guide_legend(nrow=2,byrow=TRUE)) 

g

ggsave(here::here("output/figures/Fig_stable_1.pdf"), width = 6, height = 5, dpi = 120)

Stables 2: High/Low cue by session

stable_data_df2 = stable_data  %>%
  group_by(id, trtype_str, contingency) %>%
  summarise_at(c("model_pred","prob"), mean, na.rm = TRUE)

w <- 0.4
g <- ggplot(data = stable_data_df2, aes(y = prob, x = contingency, fill = trtype_str))  
h<-c(0.6, 0.75, 0.9)
#
for (s in seq(3)) {
  g<- g + geom_segment(x=s-w,y=h[s], xend=s+w, yend=h[s] , linetype="dashed", color=pal[5], size=0.5 , alpha=0.3)
  g <- g + geom_segment(x=s-w,y=1-h[s], xend=s+w, yend=1-h[s] , linetype="dashed", color=pal[7], size=0.5 , alpha=0.3)
}
g<- g+geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, show.legend = FALSE, lwd=1.2) +
geom_point(aes(y = prob, colour=trtype_str), position = position_jitterdodge(  jitter.width = .15,  dodge.width = 0.2), size = 2, alpha = 0.8, show.legend = FALSE) +
geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7,lwd=1.2, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=4, shape=23, aes(y=model_pred, group=trtype_str), color="black", 
             position = position_dodge( 0.2), stroke=1.5, show.legend = FALSE) +
  
expand_limits(x = 3) +
scale_color_manual(values = c(pal[5], pal[7])) +
scale_fill_manual(values = c(pal[5], pal[7])) +  
theme_bw() +
raincloud_theme +
labs(y= "Mean probability", x="Session") +
guides(fill=guide_legend(nrow=2,byrow=TRUE)) 
g

ggsave(here::here("output/figures/Fig_stable_2.pdf"), width = 8, height = 5, dpi = 120)

Stables 3: High/Low cue by median-split TA

stable_data_df3 = stable_data %>%
  group_by(id, trtype_str, tabin) %>%
  summarise_at(c("model_pred","prob"), mean, na.rm = TRUE)

g <- ggplot(data = stable_data_df3, aes(y = prob, x = trtype_str, fill = interaction(tabin,trtype_str))) +
geom_hline(yintercept = 0.75, color=pal[5], linetype="longdash")  +
geom_hline(yintercept = 0.25, color=pal[7], linetype="longdash")  +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2) +

geom_point(aes(color= interaction( tabin, trtype_str)), position =   position_jitterdodge(  jitter.width = .45,  dodge.width = 0.05), size = 2, alpha = 0.5, show.legend=FALSE) + 
geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +  
stat_summary(geom="point", fun="mean", size=4, shape=23, aes(y=model_pred, group=interaction(tabin,trtype_str)), color="black", 
             position = position_dodge( 0.2), stroke=1.5, show.legend = FALSE) +  
  
expand_limits(x = 3) +
scale_color_manual(values = c(pal[6], pal[5], pal[8], pal[7])) +
scale_fill_manual(values = c(pal[6], pal[5], pal[8], pal[7]), name = "", labels = c("HC: Low Anxiety", "HC: High Anxiety", "LC: Low Anxiety", "LC: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(legend.position = "right") + 
labs(y= "Mean probability", x="Session") +
  scale_x_discrete(labels=c("hi-prob" = "HC", "low-prob" =  "LC")) + 
guides(fill=guide_legend(nrow=4,byrow=TRUE)) 
g

ggsave(here::here("output/figures/Fig_stable_3.pdf"), width = 7, height = 5, dpi = 120)
stable_data_df3 = stable_data %>%
  group_by(id, trtype_str, ta_orig_scale) %>%
  summarise_at(c("model_pred","prob"), mean, na.rm = TRUE)

g <- ggplot(data = stable_data_df3, aes(y = prob, x = ta_orig_scale, fill = trtype_str)) +
geom_hline(yintercept = 0.75, color=pal[5], linetype="longdash")  +
geom_hline(yintercept = 0.25, color=pal[7], linetype="longdash")  +
#geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2) +

geom_point(aes(color= trtype_str), position =   position_jitterdodge(  jitter.width = .45,  dodge.width = 0.05), size = 2, alpha = 0.5, show.legend=FALSE) + 
geom_smooth(aes(color= trtype_str, fill=trtype_str), method='lm', formula= y~x, alpha=0.3, show.legend = FALSE) +

#geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +  
#stat_summary(geom="point", fun="mean", size=4, shape=23, aes(y=model_pred, group=interaction(tabin,trtype_str)), color="black", 
#             position = position_dodge( 0.2), stroke=1.5, show.legend = FALSE) +  
  
expand_limits(x = 3) +
scale_color_manual(values = c(pal[5], pal[7], pal[8], pal[7])) +
scale_fill_manual(values = c(pal[5], pal[7], pal[8], pal[7]), name = "", labels = c("HC: Low Anxiety", "HC: High Anxiety", "LC: Low Anxiety", "LC: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(legend.position = "right") + 
labs(y= "Mean probability", x="Trait anxiety") +
  scale_x_discrete(labels=c("hi-prob" = "HC", "low-prob" =  "LC")) + 
guides(fill=guide_legend(nrow=4,byrow=TRUE))  + 
  xlim(20, 70)
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
g

ggsave(here::here("output/figures/Fig_stable_3_ver2.pdf"), width = 4, height = 5, dpi = 120)

Stables 4: Error from true by cue type, session and TA

stable_data_df4 = stable_data %>%
  group_by(id, trtype_str, contingency, tabin) %>%
  summarise_at(c("model_pred", "err"), mean, na.rm = TRUE)

stable_data_df4 <- within(stable_data_df4, err[trtype_str=="hi-prob"] <- (err[trtype_str=="hi-prob"]))
# https://stackoverflow.com/questions/16026215/generate-ggplot2-boxplot-with-different-colours-for-multiple-groups
g <- ggplot(data = stable_data_df4, aes(y = err, x = trtype_str, fill = interaction( tabin, trtype_str))) +
geom_hline(yintercept=c(0), linetype="dashed") +
geom_point(aes(color= interaction( tabin, trtype_str)), position =   position_jitterdodge(  jitter.width = .45,  dodge.width = 1), 
           size = 2, alpha = 0.5, show.legend=FALSE) + 
geom_boxplot(inherit.aes = TRUE, width = 0.8, lwd=1.2, outlier.shape = NA, alpha = 0.7, show.legend = FALSE) +  
stat_summary(geom="point", fun="mean", size=4, shape=23, aes(y=err, group=interaction(tabin,trtype_str)), color="black", 
             position = position_dodge( 0.8), stroke=1.5, show.legend = FALSE) +
  
expand_limits(x = 3) +
scale_color_manual(values = c(pal[6], pal[5], pal[8], pal[7])) +
scale_fill_manual(values = c(pal[6], pal[5], pal[8], pal[7]), name = "", labels = c("HC: Low Anxiety", "HC: High Anxiety", "LC: Low Anxiety", "LC: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(legend.position = "right",
      strip.text.x = element_text(size = 12,  face = "bold"),
      strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")
      ) + 
labs(y= "Error", x="Session")  +
scale_x_discrete(labels=c("hi-prob" = "HC", "low-prob" =  "LC")) + 
ylim(-0.6, 0.6) +
facet_grid(cols=vars(contingency)) 
g

ggsave(here::here("output/figures/Fig_stable_4.pdf"), width = 8, height = 5, dpi = 120)

2. Analysis of ratings - reversal cue

Stats

rev_data<-read.csv(here::here("data/reversal_cue_data.csv"))
rev_data = assign_var_types(rev_data, c("phase_str", "study_str", "id", "contingency", "ta"))


dt_count <- rev_data %>%          
  group_by(contingency) %>%
  summarise(count = n_distinct(id))
dt_count
## # A tibble: 3 × 2
##   contingency count
##   <fct>       <int>
## 1 60/40          36
## 2 75/25          88
## 3 90/10          37
m2a <-lmer(prob ~ ta*contingency*phase_str + (1|id)  + (1|study_str), rev_data) 
## boundary (singular) fit: see help('isSingular')
m2b<-lmer(prob ~ ta*contingency*phase_str + (1|id)  + (1+ta+contingency+phase_str|study_str), rev_data)
## boundary (singular) fit: see help('isSingular')
anova(m2a, m2b)
## refitting model(s) with ML (instead of REML)
## Data: rev_data
## Models:
## m2a: prob ~ ta * contingency * phase_str + (1 | id) + (1 | study_str)
## m2b: prob ~ ta * contingency * phase_str + (1 | id) + (1 + ta + contingency + phase_str | study_str)
##     npar    AIC      BIC logLik deviance  Chisq Df Pr(>Chisq)
## m2a   15 -156.7 -100.079 93.349   -186.7                     
## m2b   29 -128.7  -19.237 93.349   -186.7 0.0018 14          1
# anova
anova(m2a)
## Type III Analysis of Variance Table with Satterthwaite's method
##                          Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
## ta                       0.0454  0.0454     1  86.976   1.6244 0.2058701    
## contingency              0.0346  0.0173     2 276.994   0.6199 0.5387247    
## phase_str                4.6362  4.6362     1 234.049 165.9236 < 2.2e-16 ***
## ta:contingency           0.0453  0.0226     2 277.980   0.8103 0.4457800    
## ta:phase_str             0.3843  0.3843     1 234.049  13.7552 0.0002601 ***
## contingency:phase_str    0.4950  0.2475     2 234.049   8.8578 0.0001958 ***
## ta:contingency:phase_str 0.0922  0.0461     2 234.049   1.6495 0.1943649    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# statistical assumptions
## homogeneity of variance 
rev_data$residuals <- residuals(m2a)
rev_data$residuals_sq <- abs(residuals(m2a))^2
hist(rev_data$residuals_sq)

lev_m <- lm(residuals_sq ~ id, data=rev_data) #ANOVA of the squared residuals
anova(lev_m)
## Analysis of Variance Table
## 
## Response: residuals_sq
##            Df  Sum Sq    Mean Sq F value Pr(>F)
## id         88 0.08185 0.00093012  0.8141 0.8676
## Residuals 233 0.26620 0.00114247
## normality of residuals
ggplot(data=rev_data, aes(x=residuals)) + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

car::qqPlot(rev_data$residuals)

## [1] 284 151
# effectsize
effectsize::effectsize(anova(m2a), type="eta", alternative = "two.sided")
## # Effect Size for ANOVA (Type III)
## 
## Parameter                | Eta2 (partial) |       95% CI
## --------------------------------------------------------
## ta                       |           0.02 | [0.00, 0.11]
## contingency              |       4.46e-03 | [0.00, 0.03]
## phase_str                |           0.41 | [0.32, 0.50]
## ta:contingency           |       5.80e-03 | [0.00, 0.03]
## ta:phase_str             |           0.06 | [0.01, 0.12]
## contingency:phase_str    |           0.07 | [0.02, 0.14]
## ta:contingency:phase_str |           0.01 | [0.00, 0.05]
# post-hocs
em2a1 = emmeans(m2a, specs = pairwise ~ phase_str*contingency)
## NOTE: Results may be misleading due to involvement in interactions
em2a1$emmeans
##  phase_str contingency emmean     SE   df lower.CL upper.CL
##  acq       60/40        0.630 0.0386 7.42    0.539    0.720
##  ext       60/40        0.483 0.0386 7.42    0.393    0.573
##  acq       75/25        0.714 0.0209 4.64    0.659    0.769
##  ext       75/25        0.455 0.0209 4.64    0.400    0.510
##  acq       90/10        0.763 0.0382 7.02    0.673    0.853
##  ext       90/10        0.384 0.0382 7.02    0.293    0.474
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
em2a1$contrasts
##  contrast                  estimate     SE    df t.ratio p.value
##  (acq 60/40) - (ext 60/40)   0.1464 0.0395 228.6   3.707  0.0035
##  (acq 60/40) - (acq 75/25)  -0.0842 0.0394  80.1  -2.138  0.2790
##  (acq 60/40) - (ext 75/25)   0.1743 0.0394  80.1   4.422  0.0004
##  (acq 60/40) - (acq 90/10)  -0.1337 0.0392 229.4  -3.407  0.0100
##  (acq 60/40) - (ext 90/10)   0.2458 0.0392 229.4   6.265  <.0001
##  (ext 60/40) - (acq 75/25)  -0.2307 0.0394  80.1  -5.853  <.0001
##  (ext 60/40) - (ext 75/25)   0.0278 0.0394  80.1   0.706  0.9807
##  (ext 60/40) - (acq 90/10)  -0.2801 0.0392 229.4  -7.139  <.0001
##  (ext 60/40) - (ext 90/10)   0.0994 0.0392 229.4   2.534  0.1187
##  (acq 75/25) - (ext 75/25)   0.2585 0.0252 228.6  10.247  <.0001
##  (acq 75/25) - (acq 90/10)  -0.0495 0.0390  77.6  -1.266  0.8021
##  (acq 75/25) - (ext 90/10)   0.3301 0.0390  77.6   8.453  <.0001
##  (ext 75/25) - (acq 90/10)  -0.3079 0.0390  77.6  -7.886  <.0001
##  (ext 75/25) - (ext 90/10)   0.0716 0.0390  77.6   1.833  0.4509
##  (acq 90/10) - (ext 90/10)   0.3795 0.0389 228.6   9.759  <.0001
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 6 estimates
emstr <- summary(em2a1$contrasts)
effectsize::t_to_eta2(t=emstr$t.ratio, df_error = emstr$df, alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.06           | [0.01, 0.12]
## 0.05           | [0.00, 0.17]
## 0.20           | [0.06, 0.34]
## 0.05           | [0.01, 0.11]
## 0.15           | [0.07, 0.23]
## 0.30           | [0.14, 0.44]
## 6.19e-03       | [0.00, 0.08]
## 0.18           | [0.10, 0.27]
## 0.03           | [0.00, 0.08]
## 0.31           | [0.22, 0.40]
## 0.02           | [0.00, 0.12]
## 0.48           | [0.32, 0.60]
## 0.44           | [0.28, 0.57]
## 0.04           | [0.00, 0.16]
## 0.29           | [0.20, 0.38]
# relationship with anxiety

jt<-joint_tests(m2a, by = "phase_str")
jt 
## phase_str = acq:
##  model term     df1    df2 F.ratio p.value
##  ta               1 154.55   0.904  0.3432
##  contingency      2 112.09   5.375  0.0059
##  ta:contingency   2 229.56   0.044  0.9573
## 
## phase_str = ext:
##  model term     df1    df2 F.ratio p.value
##  ta               1 154.55   9.426  0.0025
##  contingency      2 112.09   4.760  0.0104
##  ta:contingency   2 229.56   2.319  0.1006
effectsize::F_to_eta2(f=jt$F.ratio, df=jt$df1, df_error =  jt$df2, type="eta",  alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 5.82e-03       | [0.00, 0.05]
## 0.06           | [0.01, 0.14]
## 0.09           | [0.01, 0.19]
## 0.08           | [0.00, 0.18]
## 3.83e-04       | [0.00, 0.00]
## 0.02           | [0.00, 0.06]
emtr <- emtrends(m2a, pairwise ~ phase_str, var = "ta")
## NOTE: Results may be misleading due to involvement in interactions
emtr
## $emtrends
##  phase_str ta.trend      SE  df lower.CL upper.CL
##  acq        0.00160 0.00168 155 -0.00173  0.00493
##  ext       -0.00517 0.00168 155 -0.00850 -0.00184
## 
## Results are averaged over the levels of: contingency 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast  estimate      SE  df t.ratio p.value
##  acq - ext  0.00677 0.00183 229   3.709  0.0003
## 
## Results are averaged over the levels of: contingency 
## Degrees-of-freedom method: kenward-roger
sumtr<-summary(emtr$contrasts)
sumtr
##  contrast  estimate      SE  df t.ratio p.value
##  acq - ext  0.00677 0.00183 229   3.709  0.0003
## 
## Results are averaged over the levels of: contingency 
## Degrees-of-freedom method: kenward-roger
effectsize::t_to_eta2(t=sumtr$t.ratio, df_error = sumtr$df, type="eta",  alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.06           | [0.01, 0.12]
## there is no three-way interaction but this is for visualisation purpose
emtr <- emtrends(m2a, pairwise ~ phase_str * contingency, var = "ta")
emtr
## $emtrends
##  phase_str contingency ta.trend      SE  df lower.CL upper.CL
##  acq       60/40        0.00146 0.00279 303 -0.00402  0.00694
##  ext       60/40       -0.00499 0.00279 303 -0.01047  0.00049
##  acq       75/25        0.00121 0.00196 264 -0.00264  0.00507
##  ext       75/25       -0.00193 0.00196 264 -0.00578  0.00192
##  acq       90/10        0.00213 0.00274 302 -0.00326  0.00751
##  ext       90/10       -0.00859 0.00274 302 -0.01398 -0.00320
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                   estimate      SE  df t.ratio p.value
##  (acq 60/40) - (ext 60/40)  0.006454 0.00352 229   1.835  0.4456
##  (acq 60/40) - (acq 75/25)  0.000247 0.00315 271   0.079  1.0000
##  (acq 60/40) - (ext 75/25)  0.003391 0.00315 271   1.078  0.8899
##  (acq 60/40) - (acq 90/10) -0.000664 0.00349 230  -0.190  1.0000
##  (acq 60/40) - (ext 90/10)  0.010053 0.00349 230   2.880  0.0491
##  (ext 60/40) - (acq 75/25) -0.006206 0.00315 271  -1.972  0.3610
##  (ext 60/40) - (ext 75/25) -0.003062 0.00315 271  -0.973  0.9261
##  (ext 60/40) - (acq 90/10) -0.007117 0.00349 230  -2.039  0.3236
##  (ext 60/40) - (ext 90/10)  0.003600 0.00349 230   1.031  0.9070
##  (acq 75/25) - (ext 75/25)  0.003144 0.00239 229   1.316  0.7757
##  (acq 75/25) - (acq 90/10) -0.000911 0.00310 270  -0.294  0.9997
##  (acq 75/25) - (ext 90/10)  0.009806 0.00310 270   3.159  0.0216
##  (ext 75/25) - (acq 90/10) -0.004055 0.00310 270  -1.306  0.7814
##  (ext 75/25) - (ext 90/10)  0.006662 0.00310 270   2.146  0.2668
##  (acq 90/10) - (ext 90/10)  0.010717 0.00345 229   3.103  0.0259
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 6 estimates
sumtr<-summary(emtr$contrasts)
sumtr
##  contrast                   estimate      SE  df t.ratio p.value
##  (acq 60/40) - (ext 60/40)  0.006454 0.00352 229   1.835  0.4456
##  (acq 60/40) - (acq 75/25)  0.000247 0.00315 271   0.079  1.0000
##  (acq 60/40) - (ext 75/25)  0.003391 0.00315 271   1.078  0.8899
##  (acq 60/40) - (acq 90/10) -0.000664 0.00349 230  -0.190  1.0000
##  (acq 60/40) - (ext 90/10)  0.010053 0.00349 230   2.880  0.0491
##  (ext 60/40) - (acq 75/25) -0.006206 0.00315 271  -1.972  0.3610
##  (ext 60/40) - (ext 75/25) -0.003062 0.00315 271  -0.973  0.9261
##  (ext 60/40) - (acq 90/10) -0.007117 0.00349 230  -2.039  0.3236
##  (ext 60/40) - (ext 90/10)  0.003600 0.00349 230   1.031  0.9070
##  (acq 75/25) - (ext 75/25)  0.003144 0.00239 229   1.316  0.7757
##  (acq 75/25) - (acq 90/10) -0.000911 0.00310 270  -0.294  0.9997
##  (acq 75/25) - (ext 90/10)  0.009806 0.00310 270   3.159  0.0216
##  (ext 75/25) - (acq 90/10) -0.004055 0.00310 270  -1.306  0.7814
##  (ext 75/25) - (ext 90/10)  0.006662 0.00310 270   2.146  0.2668
##  (acq 90/10) - (ext 90/10)  0.010717 0.00345 229   3.103  0.0259
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 6 estimates
effectsize::t_to_eta2(t=sumtr$t.ratio, df_error = sumtr$df, type="eta",  alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.01           | [0.00, 0.06]
## 2.29e-05       | [0.00, 0.01]
## 4.27e-03       | [0.00, 0.03]
## 1.57e-04       | [0.00, 0.01]
## 0.03           | [0.00, 0.09]
## 0.01           | [0.00, 0.05]
## 3.49e-03       | [0.00, 0.03]
## 0.02           | [0.00, 0.07]
## 4.61e-03       | [0.00, 0.04]
## 7.53e-03       | [0.00, 0.04]
## 3.19e-04       | [0.00, 0.02]
## 0.04           | [0.01, 0.09]
## 6.28e-03       | [0.00, 0.04]
## 0.02           | [0.00, 0.06]
## 0.04           | [0.01, 0.10]
# one-way t-tests against true reinforcement rates
df2 <- rev_data %>% 
    group_by(contingency, tabin, phase_str) %>%
    nest() %>% 
    mutate(tt=map(data,~t.test(.x$err))) %>%
    mutate(tidied = map(tt, tidy)) %>% 
    unnest(tidied, .drop = T)
df2$p.value.fdr = p.adjust(df2$p.value, method = "fdr")
df2$sig <- 0 
df2$sig[df2$p.value.fdr<0.05] <- 1
df2
## # A tibble: 12 × 15
## # Groups:   contingency, tabin, phase_str [12]
##    contingency phase_str tabin data  tt    estimate statistic  p.value parameter
##    <fct>       <fct>     <fct> <lis> <lis>    <dbl>     <dbl>    <dbl>     <dbl>
##  1 75/25       acq       lowA… <tib… <hte…  -0.0677    -2.79  7.50e- 3        47
##  2 75/25       ext       lowA… <tib… <hte…   0.269      9.42  2.12e-12        47
##  3 75/25       acq       hiAnx <tib… <hte…  -0.0602    -2.47  1.79e- 2        39
##  4 75/25       ext       hiAnx <tib… <hte…   0.214      6.33  1.81e- 7        39
##  5 60/40       acq       hiAnx <tib… <hte…  -0.0229    -0.662 5.17e- 1        17
##  6 60/40       ext       hiAnx <tib… <hte…   0.0936     2.12  4.91e- 2        17
##  7 90/10       acq       hiAnx <tib… <hte…  -0.103     -2.13  4.80e- 2        17
##  8 90/10       ext       hiAnx <tib… <hte…   0.194      3.38  3.59e- 3        17
##  9 60/40       acq       lowA… <tib… <hte…  -0.0212    -0.502 6.22e- 1        17
## 10 60/40       ext       lowA… <tib… <hte…   0.232      5.25  6.47e- 5        17
## 11 90/10       acq       lowA… <tib… <hte…  -0.125     -3.01  7.47e- 3        18
## 12 90/10       ext       lowA… <tib… <hte…   0.351      5.99  1.16e- 5        18
## # … with 6 more variables: conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.value.fdr <dbl>, sig <dbl>
# binarized anxiety to get means per group for descriptive tables (in supplementary)
m2c <-lmer(prob ~ tabin*contingency*phase_str + (1|id)  + (1|study_str), rev_data) 
## boundary (singular) fit: see help('isSingular')
em2c1 = emmeans(m2c, specs = pairwise ~ phase_str*tabin)
## NOTE: Results may be misleading due to involvement in interactions
em2c1$emmeans
##  phase_str tabin  emmean     SE   df lower.CL upper.CL
##  acq       lowAnx  0.698 0.0303 6.03    0.624    0.772
##  ext       lowAnx  0.487 0.0303 6.03    0.412    0.561
##  acq       hiAnx   0.707 0.0322 6.00    0.628    0.786
##  ext       hiAnx   0.388 0.0322 6.00    0.309    0.467
## 
## Results are averaged over the levels of: contingency 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
em2c2 = emmeans(m2c, specs = pairwise ~ phase_str*contingency*tabin)
em2c2$emmeans
##  phase_str contingency tabin  emmean     SE   df lower.CL upper.CL
##  acq       60/40       lowAnx  0.629 0.0499 28.8    0.527    0.731
##  ext       60/40       lowAnx  0.536 0.0499 28.8    0.434    0.638
##  acq       75/25       lowAnx  0.712 0.0279 16.6    0.653    0.771
##  ext       75/25       lowAnx  0.469 0.0279 16.6    0.410    0.528
##  acq       90/10       lowAnx  0.753 0.0489 26.3    0.653    0.854
##  ext       90/10       lowAnx  0.455 0.0489 26.3    0.354    0.555
##  acq       60/40       hiAnx   0.632 0.0504 27.2    0.528    0.735
##  ext       60/40       hiAnx   0.422 0.0504 27.2    0.318    0.525
##  acq       75/25       hiAnx   0.715 0.0314 19.3    0.649    0.780
##  ext       75/25       hiAnx   0.441 0.0314 19.3    0.376    0.507
##  acq       90/10       hiAnx   0.774 0.0504 27.2    0.671    0.878
##  ext       90/10       hiAnx   0.301 0.0504 27.2    0.197    0.404
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
em2c3 = emmeans(m2c, specs = pairwise ~ phase_str*contingency)
## NOTE: Results may be misleading due to involvement in interactions
em2c3$emmeans
##  phase_str contingency emmean     SE   df lower.CL upper.CL
##  acq       60/40        0.630 0.0391 7.34    0.539    0.722
##  ext       60/40        0.479 0.0391 7.34    0.387    0.571
##  acq       75/25        0.713 0.0212 4.80    0.658    0.769
##  ext       75/25        0.455 0.0212 4.80    0.400    0.510
##  acq       90/10        0.764 0.0388 7.00    0.672    0.856
##  ext       90/10        0.378 0.0388 7.00    0.286    0.469
## 
## Results are averaged over the levels of: tabin 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Plots

reversal-locked ratings by state and anxiety

rev_data_cont <- read.csv(here::here("data/reversal_aligned_data.csv"))
rev_data_cont = assign_var_types(rev_data_cont, c("phase_str", "study_str", "id", "contingency", "ta", "half_str"))

#first summarize by participant
rev_data_cont_df = rev_data_cont %>%
  group_by(phase_str, contingency, id, trno) %>%
  summarise_at("prob", funs(mean),na.rm = TRUE)

df = rev_data_cont_df %>%
  group_by(phase_str, trno) %>%
  summarise_at("prob", funs(mean,std.error),na.rm = TRUE)
df$lower = df$mean - df$std.error
df$upper = df$mean + df$std.error

g <- ggplot(data = df, aes(y = mean, x = trno, fill=phase_str)) +
  geom_segment(x=0.5,y=0.75, xend=15, yend=0.75 , linetype="dashed", color=pal[1], size=0.5 , alpha=0.9) +
  geom_segment(x=0.5,y=0.25, xend=15, yend=0.25 , linetype="dashed", color=pal[3], size=0.5 , alpha=0.3) +
  geom_segment(x=-5,y=0.75, xend=-0.5, yend=0.75 , linetype="dashed", color=pal[3], size=0.5 , alpha=0.9) +
  geom_segment(x=-5,y=0.25, xend=-0.5, yend=0.25 , linetype="dashed", color=pal[1], size=0.5 , alpha=0.3) +
  
  geom_vline(xintercept=0, linetype="dashed") +
  geom_line(size=1.5, show.legend = TRUE) +
  geom_ribbon(aes(ymin=lower, ymax=upper),  na.rm = TRUE,alpha=0.8, show.legend = TRUE) + 
  geom_rect(xmin=10,xmax=15,ymin=0,ymax=1, size=2, color=pal2[7], fill=rgb(0.9, 0.9, 0.9), alpha=0.01) +
  scale_color_manual(values = c(pal[1], pal[3])) +
  scale_fill_manual(values = c(pal[1], pal[3]), name = "Switch type", labels = c("Low-to-high (L2H)", "High-to-low (H2L)")) +
  scale_x_continuous(breaks = c(-5, -1, 1, 5, 10, 15)) +
  ylim(0,1) + 
  theme_bw() +
  labs(y= "Rating", x="Trial (locked to reversal)")  +
  raincloud_theme  +
  theme(legend.position = c(0.35,0.9))
g

ggsave(here::here("output/figures/Fig_reversal_1.pdf"), width = 5, height = 5, dpi = 120)

reversal-locked ratings by state and anxiety

#first summarize by participant
rev_data_cont_df2 = rev_data_cont %>%
  group_by(phase_str, tabin, contingency, id, trno) %>%
  summarise_at("prob", funs(mean),na.rm = TRUE)

df = rev_data_cont_df2 %>%
  group_by(phase_str,tabin, trno) %>%
  summarise_at("prob", funs(mean,std.error),na.rm = TRUE)
df$lower = df$mean - df$std.error
df$upper = df$mean + df$std.error


g <- ggplot(data = df, aes(y = mean, x = trno, fill=interaction(phase_str, tabin))) +
  geom_segment(x=0.5,y=0.75, xend=15, yend=0.75 , linetype="dashed", color=pal[1], size=0.5 , alpha=0.9) +
  geom_segment(x=0.5,y=0.25, xend=15, yend=0.25 , linetype="dashed", color=pal[3], size=0.5 , alpha=0.3) +
  geom_segment(x=-5,y=0.75, xend=-0.5, yend=0.75 , linetype="dashed", color=pal[3], size=0.5 , alpha=0.9) +
  geom_segment(x=-5,y=0.25, xend=-0.5, yend=0.25 , linetype="dashed", color=pal[1], size=0.5 , alpha=0.3) +
  geom_vline(xintercept=0, linetype="dashed") +
  geom_line(size=1.5, show.legend = TRUE) +
  geom_ribbon(aes(ymin=lower, ymax=upper),  na.rm = TRUE,alpha=0.6,show.legend = TRUE) + 
  scale_fill_manual(values = c(pal[2], pal[4], pal[1], pal[3]), name = "", labels = c("Low Anx: L2H", "Low Anx: H2L", "High Anx: L2H", "High Anx: H2L")) +
  scale_color_manual(values = c(pal[2], pal[4], pal[1], pal[3])) + 
  ylim(0,1) + 
  scale_x_continuous(breaks = c(-5, -1, 1, 5, 10), limits=c(-5,10)) +
  theme_bw() +
  labs(y= "Rating", x="Trial (locked to reversal)")  +
  raincloud_theme  +
  theme(legend.position = c(0.75,0.25))  
g

ggsave(here::here("output/figures/Fig_reversal_2.pdf"), width = 5, height = 5, dpi = 120)

reversal cue by anxiety, session and a state

df = rev_data %>%
  group_by(id, phase_str, contingency, tabin) %>%
  summarise_at("prob", mean, na.rm = TRUE)

df_levs <- data.frame(contingency = c("60/40", "60/40", "75/25", "75/25", "90/10", "90/10"), 
                      phase_str = c("acq", "ext", "acq", "ext", "acq", "ext"),
                      level = c(60, 40, 75, 25, 90,10))
g <- ggplot(data = df, aes(y = prob, x = phase_str, fill = interaction( tabin, phase_str))) +
  geom_hline(data = df_levs[df_levs$phase_str == "acq",], aes(yintercept = level/100), color=pal[1], linetype="longdash", alpha =0.8) +
  geom_hline(data = df_levs[df_levs$phase_str == "ext",], aes(yintercept = level/100), color=pal[3], linetype="longdash", alpha =0.8) +
geom_point(aes(color= interaction( tabin, phase_str)), position =   position_jitterdodge(  jitter.width = .45,  dodge.width = 1), size = 2, alpha = 0.5, show.legend=FALSE) + 
geom_boxplot(inherit.aes = TRUE, width = 0.8, lwd=1.2, outlier.shape = NA, alpha = 0.7) + 
stat_summary(geom="point", fun="mean", size=4, shape=23, aes(group=interaction(tabin,phase_str)), color="black", 
             position = position_dodge( 0.8), stroke=1.5, show.legend = FALSE) +  
  
expand_limits(x = 3) +
scale_color_manual(values = c(pal[2], pal[1], pal[4], pal[3])) +
scale_fill_manual(values = c(pal[2], pal[1], pal[4], pal[3]), name = "", labels = c("HS: Low Anxiety", "HS: High Anxiety", "LS: Low Anxiety", "LS: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14,  face = "bold"),
      strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
theme(legend.position = "top") +  #c(0.8, 0.2)
labs(y= "Mean rating", x="State")  +
scale_x_discrete(labels=c("acq" = "High", "ext" =  "Low")) + 
facet_grid(cols=vars(contingency)) +
  guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
ylim(0,1)
g

ggsave(here::here("output/figures/Fig_reversal_3.pdf"), width = 6, height = 5, dpi = 120)

without session to better reflect the statistical model

df = rev_data %>%
  group_by(id, phase_str, contingency, ta_orig_scale) %>%
  summarise_at("prob", mean, na.rm = TRUE)

df_levs <- data.frame(contingency = c("75/25", "75/25"), 
                      phase_str = c("acq", "ext", "acq", "ext", "acq", "ext"),
                      level = c(75, 25))
g <- ggplot(data = df, aes(y = prob, x = ta_orig_scale, fill = interaction(phase_str))) +
geom_hline(data = df_levs[df_levs$phase_str == "acq",], aes(yintercept = level/100), color=pal[1], linetype="longdash", alpha =0.8) +
geom_hline(data = df_levs[df_levs$phase_str == "ext",], aes(yintercept = level/100), color=pal[3], linetype="longdash", alpha =0.8) +
geom_point(aes(color= phase_str), position =   position_jitterdodge(  jitter.width = .45,  dodge.width = 1), size = 2, alpha = 0.5, show.legend=FALSE) + 
geom_smooth(aes(color= phase_str, fill=phase_str), method='lm', formula= y~x, alpha=0.3, show.legend = FALSE) +
#geom_boxplot(inherit.aes = TRUE, width = 0.8, lwd=1.2, outlier.shape = NA, alpha = 0.7) + 
#stat_summary(geom="point", fun="mean", size=4, shape=23, aes(group=interaction(tabin,phase_str)), color="black", 
#             position = position_dodge( 0.8), stroke=1.5, show.legend = FALSE) +  
  
expand_limits(x = 3) +
scale_color_manual(values = c(pal[1], pal[3], pal[4], pal[3])) +
scale_fill_manual(values = c(pal[1], pal[3], pal[4], pal[3]), name = "", labels = c("HS: Low Anxiety", "HS: High Anxiety", "LS: Low Anxiety", "LS: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14,  face = "bold"),
      strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
theme(legend.position = "top") +  #c(0.8, 0.2)
labs(y= "Mean rating on trials 10+", x="Trait anxiety")  +
#scale_x_discrete(labels=c("acq" = "High", "ext" =  "Low")) + 
#facet_grid(cols=vars(contingency)) +
  guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
ylim(0,1) + 
  xlim(20,70)
g

ggsave(here::here("output/figures/Fig_reversal_3_collapsed.pdf"), width = 4, height = 5, dpi = 120)

Only 75/25 by study

df = rev_data %>%
  group_by(id, phase_str, contingency, tabin,ta, study_str) %>%
  summarise_at("prob", mean, na.rm = TRUE)
df <- df[df$contingency %in% "75/25",]

m<-lmer(prob ~ phase_str*study_str + (1|id) , df)
anova(m)
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF DenDF  F value Pr(>F)    
## phase_str           2.87820 2.87820     1    85 116.2776 <2e-16 ***
## study_str           0.03662 0.01831     2    85   0.7398 0.4803    
## phase_str:study_str 0.05418 0.02709     2    85   1.0945 0.3394    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
joint_tests(m, by = "phase_str")
## phase_str = acq:
##  model term df1    df2 F.ratio p.value
##  study_str    2 158.14   0.351  0.7047
## 
## phase_str = ext:
##  model term df1    df2 F.ratio p.value
##  study_str    2 158.14   1.386  0.2530
df_levs <- data.frame(contingency = c("75/25", "75/25", "75/25", "75/25", "75/25", "75/25"), 
                      phase_str = c("acq", "ext", "acq", "ext", "acq", "ext"),
                      level = c(75, 25, 75, 25, 75, 25))
g <- ggplot(data = df, aes(y = prob, x = phase_str, fill = interaction( tabin, phase_str))) +
  geom_hline(data = df_levs[df_levs$phase_str == "acq",], aes(yintercept = level/100), color=pal[1], linetype="longdash", alpha =0.8) +
  geom_hline(data = df_levs[df_levs$phase_str == "ext",], aes(yintercept = level/100), color=pal[3], linetype="longdash", alpha =0.8) +
geom_point(aes(color= interaction( tabin, phase_str)), position =   position_jitterdodge(  jitter.width = .45,  dodge.width = 1), size = 2, alpha = 0.5, show.legend=FALSE) + 
geom_boxplot(inherit.aes = TRUE, width = 0.8, lwd=1.2, outlier.shape = NA, alpha = 0.7) + 

  
expand_limits(x = 3) +
scale_color_manual(values = c(pal[2], pal[1], pal[4], pal[3])) +
scale_fill_manual(values = c(pal[2], pal[1], pal[4], pal[3]), name = "", labels = c("HS: Low Anxiety", "HS: High Anxiety", "LS: Low Anxiety", "LS: High Anxiety")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14,  face = "bold"),
      strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
theme(legend.position = "top") +  #c(0.8, 0.2)
labs(y= "Mean rating", x="State")  +
scale_x_discrete(labels=c("acq" = "High", "ext" =  "Low")) + 
facet_grid(cols=vars(study_str)) +
  guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
ylim(0,1)
g

ggsave(here::here("output/figures/SuppFig5_raw.pdf"), width = 6, height = 5, dpi = 120)

3. Analysis of slopes and switchpoints

Stats

Slopes by anxiety and state

rev_data <- read.csv(here::here("data/full_REV_data.csv"))

# We'll be fitting the slopes to the first ten trials after reversal
rev_data<-rev_data[which(rev_data$trno %in% seq(10)),]
rev_data = assign_var_types(rev_data, c("phase_str", "trno", "study_str", "id", "contingency", "ta", "half_str"))


# prepare data set
conts <- c("60/40", "75/25", "90/10")
phases <- c("acq", "ext") # acq= high-state; ext = low-state
halves <- c("first", "second")
df <- rev_data[,c("id", "ta")]
df<-df[!duplicated(df$id),]
df["tabin"]<-dicho(df$ta)
df$tabin <- mapvalues(df$tabin,
                           from = c(0,1),
                           to = c("lowAnx", "hiAnx"))
df["tabin"]<-to_factor(df["tabin"])

# run lmer for each combination of session, phase and half
dt <- data.frame()
for (c in conts) {
  for (ph in phases) {
    for (h in halves)  {
      adata <- rev_data[which((rev_data$phase_str %in% ph) & (rev_data$contingency %in% c) & (rev_data$half_str %in% h)),]
      mm = lmer(prob ~ 1 + (trno|id) , adata) # This estimates slope for each participant
      rn_a<- data.frame(ranef(mm))
      rn_a["id"] <- rn_a$grp
      df_a<-data.frame()
      df_a <- merge(df, rn_a[rn_a$term=="trno",], by="id")
      df_a["phase"] <- ph
      df_a["contingency"] <- c
      df_a["half"] <- h
      dt <- rbind(dt, df_a)
    }
  }
}
dt$phase_str = to_factor(dt$phase)
dt$contingency = to_factor(dt$contingency)
dt$half_str = to_factor(dt$half)
dt$abscondval <- abs(dt$condval)

# Also save the slope estiamtes for other anlayses
dt$slope <- dt$abscondval
write.csv(dt, here::here("data", "slope_estimates.csv"))

# Ns
N = dt %>%
  group_by(contingency) %>%
  count()

dt_count <- dt %>%          
  group_by(contingency) %>%
  summarise(count = n_distinct(id))
dt_count
## # A tibble: 3 × 2
##   contingency count
##   <fct>       <int>
## 1 60/40          36
## 2 75/25          88
## 3 90/10          37
# Test on slopes (absolute slopes)
dfst<-read.csv(here::here("data/model_fit_final_full.csv"))
dfst<-dfst[,c("ids", "study", "contingency")]
dfst<-assign_var_types(dfst, c( "study", "ids")) 
dfst$id <- dfst$ids
dfst<-dfst[,c("study_str", "id", "contingency")]
dt <- base::merge(dt, dfst, by=c("id", "contingency"), all.x=TRUE)



# Report mean slopes using non-absolute values
m3b<-lmer(condval ~ ta*contingency*phase_str + (1|id) + (1|study_str), dt)
## boundary (singular) fit: see help('isSingular')
em3b = emmeans(m3b, specs = pairwise ~ phase_str)
## NOTE: Results may be misleading due to involvement in interactions
em3b$emmeans
##  phase_str  emmean      SE   df lower.CL upper.CL
##  acq        0.0244 0.00287 0.75  -0.0578   0.1067
##  ext       -0.0233 0.00290 0.74  -0.1098   0.0632
## 
## Results are averaged over the levels of: contingency 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
jt<-joint_tests(m3b, by = "contingency")
jt
## contingency = 60/40:
##  model term   df1    df2 F.ratio p.value
##  ta             1 290.73   0.901  0.3434
##  phase_str      1 554.55  33.404  <.0001
##  ta:phase_str   1 554.55   4.629  0.0319
## 
## contingency = 75/25:
##  model term   df1    df2 F.ratio p.value
##  ta             1 292.03   0.518  0.4725
##  phase_str      1 561.13 141.254  <.0001
##  ta:phase_str   1 562.70   0.050  0.8234
## 
## contingency = 90/10:
##  model term   df1    df2 F.ratio p.value
##  ta             1 290.73   0.005  0.9441
##  phase_str      1 554.55 316.719  <.0001
##  ta:phase_str   1 554.55  29.456  <.0001
# Main statistical model using absolute value of the slope
m3a<-lmer(abscondval ~ ta*contingency*phase_str + (1|id) + (1|study_str), dt)
anova(m3a)
## Type III Analysis of Variance Table with Satterthwaite's method
##                            Sum Sq   Mean Sq NumDF  DenDF F value    Pr(>F)    
## ta                       0.002629 0.0026291     1  87.42  5.8596  0.017561 *  
## contingency              0.048342 0.0241710     2  96.11 53.8716 < 2.2e-16 ***
## phase_str                0.000531 0.0005307     1 541.05  1.1828  0.277262    
## ta:contingency           0.004787 0.0023936     2 599.71  5.3348  0.005052 ** 
## ta:phase_str             0.000161 0.0001609     1 542.26  0.3587  0.549489    
## contingency:phase_str    0.000383 0.0001914     2 541.29  0.4267  0.652896    
## ta:contingency:phase_str 0.000856 0.0004278     2 542.73  0.9535  0.386053    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# statistical assumptions
## homogeneity of variance 
dt$residuals <- residuals(m3a)
dt$residuals_sq <- abs(residuals(m3a))^2
hist(dt$residuals_sq)

lev_m <- lm(residuals_sq ~ id, data=dt) #ANOVA of the squared residuals
anova(lev_m)
## Analysis of Variance Table
## 
## Response: residuals_sq
##            Df     Sum Sq    Mean Sq F value Pr(>F)
## id         88 5.0268e-05 5.7123e-07  1.2109  0.107
## Residuals 544 2.5662e-04 4.7172e-07
## normality of residuals
ggplot(data=dt, aes(x=residuals)) + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

car::qqPlot(dt$residuals)

## [1] 244 593
# effectsize
effectsize::effectsize(anova(m3a), type="eta", alternative = "two.sided")
## # Effect Size for ANOVA (Type III)
## 
## Parameter                | Eta2 (partial) |       95% CI
## --------------------------------------------------------
## ta                       |           0.06 | [0.00, 0.18]
## contingency              |           0.53 | [0.39, 0.63]
## phase_str                |       2.18e-03 | [0.00, 0.02]
## ta:contingency           |           0.02 | [0.00, 0.04]
## ta:phase_str             |       6.61e-04 | [0.00, 0.01]
## contingency:phase_str    |       1.57e-03 | [0.00, 0.01]
## ta:contingency:phase_str |       3.50e-03 | [0.00, 0.02]
# # Post-hocs
em3a = emmeans(m3a, specs = pairwise ~ contingency)
## NOTE: Results may be misleading due to involvement in interactions
em3a$emmeans
##  contingency emmean      SE   df lower.CL upper.CL
##  60/40       0.0200 0.00305 3.72   0.0113   0.0287
##  75/25       0.0252 0.00188 1.85   0.0165   0.0340
##  90/10       0.0445 0.00303 3.57   0.0357   0.0534
## 
## Results are averaged over the levels of: phase_str 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
em3a$contrasts
##  contrast          estimate      SE    df t.ratio p.value
##  (60/40) - (75/25) -0.00521 0.00271  88.2  -1.921  0.1387
##  (60/40) - (90/10) -0.02452 0.00249 541.2  -9.844  <.0001
##  (75/25) - (90/10) -0.01931 0.00269  84.9  -7.191  <.0001
## 
## Results are averaged over the levels of: phase_str 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
emstr <- summary(em3a$contrasts)
effectsize::t_to_eta2(t=emstr$t.ratio, df_error = emstr$df, alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.04           | [0.00, 0.15]
## 0.15           | [0.10, 0.21]
## 0.38           | [0.22, 0.51]
# # Interactions with trait anxiety
trdf<-emtrends(m3a, pairwise ~ contingency, var = "ta")
## NOTE: Results may be misleading due to involvement in interactions
trdf
## $emtrends
##  contingency ta.trend       SE  df  lower.CL upper.CL
##  60/40       0.000185 0.000208 241 -0.000224 0.000593
##  75/25       0.000110 0.000156 141 -0.000198 0.000417
##  90/10       0.000748 0.000205 233  0.000345 0.001151
## 
## Results are averaged over the levels of: phase_str 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate       SE  df t.ratio p.value
##  (60/40) - (75/25)  0.000075 0.000213 620   0.353  0.9336
##  (60/40) - (90/10) -0.000564 0.000222 542  -2.544  0.0302
##  (75/25) - (90/10) -0.000639 0.000209 620  -3.049  0.0067
## 
## Results are averaged over the levels of: phase_str 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
emstr <- summary(trdf$contrasts)
effectsize::t_to_eta2(t=emstr$t.ratio, df_error = emstr$df, alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 2.01e-04       | [0.00, 0.01]
## 0.01           | [0.00, 0.04]
## 0.01           | [0.00, 0.04]
jt<-joint_tests(m3a, by = "contingency")
jt
## contingency = 60/40:
##  model term   df1    df2 F.ratio p.value
##  ta             1 240.53   0.791  0.3746
##  phase_str      1 538.11   0.040  0.8416
##  ta:phase_str   1 538.11   0.453  0.5011
## 
## contingency = 75/25:
##  model term   df1    df2 F.ratio p.value
##  ta             1 141.34   0.496  0.4823
##  phase_str      1 545.42   0.030  0.8622
##  ta:phase_str   1 547.94   2.041  0.1536
## 
## contingency = 90/10:
##  model term   df1    df2 F.ratio p.value
##  ta             1 232.86  13.385  0.0003
##  phase_str      1 538.11   0.691  0.4063
##  ta:phase_str   1 538.11   0.379  0.5385
jt$F.ratio
## [1]  0.791  0.496 13.385  0.040  0.030  0.691  0.453  2.041  0.379
effectsize::F_to_eta2(f=jt$F.ratio , df=jt$df1, df_error = jt$df2, alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 3.28e-03       | [0.00, 0.03]
## 3.50e-03       | [0.00, 0.05]
## 0.05           | [0.01, 0.12]
## 7.43e-05       | [0.00, 0.01]
## 5.50e-05       | [0.00, 0.01]
## 1.28e-03       | [0.00, 0.01]
## 8.41e-04       | [0.00, 0.01]
## 3.71e-03       | [0.00, 0.02]
## 7.04e-04       | [0.00, 0.01]
jt2<-joint_tests(m3a, by = c("contingency", "phase_str"))
jt2
## contingency = 60/40, phase_str = acq:
##  model term df1    df2 F.ratio p.value
##  ta           1 435.34   1.244  0.2652
## 
## contingency = 75/25, phase_str = acq:
##  model term df1    df2 F.ratio p.value
##  ta           1 283.04   0.068  0.7950
## 
## contingency = 90/10, phase_str = acq:
##  model term df1    df2 F.ratio p.value
##  ta           1 425.79   6.484  0.0112
## 
## contingency = 60/40, phase_str = ext:
##  model term df1    df2 F.ratio p.value
##  ta           1 435.34   0.091  0.7634
## 
## contingency = 75/25, phase_str = ext:
##  model term df1    df2 F.ratio p.value
##  ta           1 292.30   1.946  0.1641
## 
## contingency = 90/10, phase_str = ext:
##  model term df1    df2 F.ratio p.value
##  ta           1 425.79  10.817  0.0011

slope and relative model fit

fit_data <- read.csv(here::here("data/model_fit_final_full.csv"))
fit_data = assign_var_types(fit_data, c("contingency", "ta", "study", "ids"))
fit_data$id <- fit_data$ids
dt$slope <- abs(dt$condval)
dt2 = dt %>%
  group_by(id, contingency) %>%
  summarise_at(c("slope"), mean, na.rm = TRUE)

write.csv(dt2, here::here("../output", "tables_and_csvs", "est_slopes_1-10.csv"), row.names = FALSE)

df <- join(x=as.data.frame(dt2[,c("contingency", "slope", "id")]), y=as.data.frame(fit_data[, c("contingency", "id", "m1m3_diff", "ta")]), by=c("contingency", "id"))

### correlation between slope and relative model fit
df = df %>%
  group_by(id ) %>%
  summarise_at(c("slope", "m1m3_diff"), mean, na.rm = TRUE)
df<- as.data.frame(df)
res<-cor.test(df$slope, df$m1m3_diff)
res
## 
##  Pearson's product-moment correlation
## 
## data:  df$slope and df$m1m3_diff
## t = 3.6114, df = 87, p-value = 0.0005087
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1652314 0.5295049
## sample estimates:
##       cor 
## 0.3610637
effectsize::t_to_eta2(t=res$statistic, df_error = res$parameter)
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.13           | [0.04, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Plots

Slopes by anxiety and state

g <- ggplot(data = dt, aes(y = condval, x = phase, fill = interaction(tabin, phase_str))) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .7,lwd=1.2) +
geom_point(aes(color= interaction( tabin, phase_str)), position = position_jitterdodge(  jitter.width = .15,  dodge.width = 0.2), size = 2, alpha = 0.5, show.legend=FALSE) + geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +
geom_hline(yintercept=c(0), linetype="dashed") +
stat_summary(geom="point", fun="mean", size=3, shape=23, aes(group=interaction(tabin,phase_str)), color="black", 
             position = position_dodge( 0.2), stroke=1, show.legend = FALSE) +    
  
  
expand_limits(x = 3) +
scale_color_manual(values = c(pal[2], pal[1], pal[4], pal[3])) +
scale_fill_manual(values = c(pal[2], pal[1], pal[4], pal[3]), name = "Trait anxiety", labels = c("Low", "High")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14,  face = "bold"),
    strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
labs(y= "Slope on trials 1 - 10", x="State") + 
scale_x_discrete(labels=c("acq" = "High", "ext" =  "Low")) +
scale_y_continuous(labels = function(x) paste0(x*100, "%")) +
facet_grid(cols=vars(contingency)) 
g

ggsave(here::here("output/figures/Fig_slopes_1.pdf"), width = 10, height = 3, dpi = 120)

Correlation of TA and slope in 90/10 session

### assumes previous dt
dt$slope <- abs(dt$condval)
dt2 = dt %>%
  group_by(id, contingency) %>%
  summarise_at(c("slope"), mean, na.rm = TRUE)

df <- join(x=as.data.frame(dt2[,c("contingency", "slope", "id")]), y=as.data.frame(fit_data[, c("contingency", "id", "m1m3_diff", "ta_orig_scale")]), by=c("contingency", "id"))

df<-df[df$contingency %in% "90/10",]
g <- ggplot(data = df, aes(x = ta_orig_scale, y = slope ), show.legend = FALSE) +
geom_point( ) +
geom_smooth(method='lm', formula= y~x, alpha=0.3, show.legend = FALSE, color="black") +
scale_color_gradient(low=pal[9], high=pal[9]) +
theme_bw() +
raincloud_theme +
  stat_cor(method = "pearson",  alternative = "two.sided",  cor.coef.name = c("r"),  label.sep = ", ",  label.x.npc = "left",  label.y.npc = "top", digits = 3, r.digits = 3, p.digits = 3) + 
  labs(x= "Trait anxiety", y="Slope on trials 1-10")   
g

ggsave(here::here("output/figures/Fig_slopes_2.pdf"), width = 5, height = 4, dpi = 120)

Switch steepness and switchpoint

stsw_data <- read.csv(here::here("data/steepness_and_switchpoint_full.csv"))
stsw_data = assign_var_types(stsw_data, c("study", "ta", "phase_str", "half_str", "cont", "tabin", "session"))

stsw_data = stsw_data %>%
  group_by(id, cont,phase_str, ta, study_str) %>%
  summarise_at(c("log_steepness", "steepness","switchpoint"), mean, na.rm = TRUE)

m <- lmer(log_steepness ~ ta*cont*phase_str + (1|id) + (1|study_str) , stsw_data)
anova(m)
## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## ta                 5.6574  5.6574     1  94.376  7.2956  0.008195 ** 
## cont              31.0785 15.5392     2 205.307 20.0387 1.122e-08 ***
## phase_str          2.1806  2.1806     1 226.858  2.8120  0.094941 .  
## ta:cont            2.3203  1.1602     2 260.232  1.4961  0.225924    
## ta:phase_str       0.1585  0.1585     1 226.858  0.2045  0.651579    
## cont:phase_str     1.1010  0.5505     2 226.858  0.7099  0.492788    
## ta:cont:phase_str  0.1640  0.0820     2 226.858  0.1057  0.899727    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effectsize
effectsize::effectsize(anova(m), type="eta")
## # Effect Size for ANOVA (Type III)
## 
## Parameter         | Eta2 (partial) |       95% CI
## -------------------------------------------------
## ta                |           0.07 | [0.01, 1.00]
## cont              |           0.16 | [0.09, 1.00]
## phase_str         |           0.01 | [0.00, 1.00]
## ta:cont           |           0.01 | [0.00, 1.00]
## ta:phase_str      |       9.00e-04 | [0.00, 1.00]
## cont:phase_str    |       6.22e-03 | [0.00, 1.00]
## ta:cont:phase_str |       9.31e-04 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].
col<- "slategray"
### Steepness plot
g <- ggplot(aes(x=cont, y=log_steepness, fill=cont), data=stsw_data) + 
  geom_flat_violin(position = position_nudge(x = .2, y = 0), lwd=1.2, alpha = .8, show.legend = FALSE) +
  geom_point(aes(y = log_steepness, color=cont), position = position_jitter(width = .15), size = 2, alpha = 0.5, show.legend = FALSE) +
  geom_boxplot(width = .2, outlier.shape = NA, lwd=1.2,alpha = 0.7, show.legend = FALSE) + 
  stat_summary(geom="point", fun="mean", size=5, shape=23, aes(group=cont), color="black", 
             position = position_dodge( 0), stroke=1, show.legend = FALSE) + 
  expand_limits(x = 3) +
  scale_color_manual(values=c(col, col, col)) +
  scale_fill_manual(values=c(col, col, col)) +
  # coord_flip() +
  theme_bw() +
  raincloud_theme +
    labs(y= "Estimated (log) Steepness", x="Session") 
  g

ggsave(here::here("output/figures/Fig_steepness.pdf"), width = 5, height = 4, dpi = 120)


### Switchpoitn plot
col <- "slategray"
g <- ggplot(aes(x=cont, y=switchpoint, fill=cont), data=stsw_data) + 
  geom_flat_violin(position = position_nudge(x = .2, y = 0), lwd=1.2, alpha = .8, show.legend = FALSE) +
  geom_point(aes(y = switchpoint, color=cont), position = position_jitter(width = .15), size = 2, alpha = 0.5, show.legend = FALSE) +
  geom_boxplot(width = .2, outlier.shape = NA, lwd=1.2, alpha = 0.7, show.legend = FALSE) +
  stat_summary(geom="point", fun="mean", size=5, shape=23, aes(group=cont), color="black", 
             position = position_dodge( 0), stroke=1, show.legend = FALSE) + 
  
  expand_limits(x = 3) +
  scale_color_manual(values=c(col, col, col)) +
  scale_fill_manual(values=c(col, col, col)) +
  theme_bw() +
  raincloud_theme +
  labs(y= "Estimated switch point", x="Session") 
  g

ggsave(here::here("output/figures/Fig_switchpoint.pdf"), width = 5, height = 4, dpi = 120)

4. Analysis of model fit and anxiety

Stats

fit_data <- read.csv(here::here("data/model_fit_final_full.csv"))
fit_data = assign_var_types(fit_data, c("study", "ta","contingency", "ids" ))

# m1m3_diff is relative model fit
m4a<-lmer(data=fit_data, m1m3_diff ~ ta*contingency + (1|ids) + (1|study_str))
## boundary (singular) fit: see help('isSingular')
anova(m4a)
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF   DenDF F value   Pr(>F)   
## ta              659.7  659.72     1  80.511  2.1302 0.148314   
## contingency    1413.6  706.82     2 104.817  2.2823 0.107102   
## ta:contingency 3217.3 1608.65     2 104.876  5.1942 0.007064 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# statistical assumptions
## homogeneity of variance 
fit_data$residuals <- residuals(m4a)
fit_data$residuals_sq <- abs(residuals(m4a))^2
lev_m <- lm(residuals_sq ~ ids, data=fit_data) #ANOVA of the squared residuals
anova(lev_m)
## Analysis of Variance Table
## 
## Response: residuals_sq
##           Df   Sum Sq Mean Sq F value Pr(>F)
## ids       88 20178129  229297  1.2707  0.147
## Residuals 72 12992832  180456
## normality of residuals
ggplot(data=fit_data, aes(x=residuals)) + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

car::qqPlot(fit_data$residuals)

## [1] 89 68
# effectsize
effectsize::effectsize(anova(m4a), type="eta", alternative = "two.sided")
## # Effect Size for ANOVA (Type III)
## 
## Parameter      | Eta2 (partial) |       95% CI
## ----------------------------------------------
## ta             |           0.03 | [0.00, 0.13]
## contingency    |           0.04 | [0.00, 0.13]
## ta:contingency |           0.09 | [0.01, 0.20]
# fit*correlation by session (p-vals uncorrected!)
jt<-joint_tests(m4a, by = "contingency")
jt 
## contingency = 60/40:
##  model term df1    df2 F.ratio p.value
##  ta           1 153.15   0.027  0.8691
## 
## contingency = 75/25:
##  model term df1    df2 F.ratio p.value
##  ta           1 143.10   0.098  0.7542
## 
## contingency = 90/10:
##  model term df1    df2 F.ratio p.value
##  ta           1 153.06   9.611  0.0023
effectsize::F_to_eta2(f=jt$F.ratio, df=jt$df1, df_error =  jt$df2, type="eta",  alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 1.76e-04       | [0.00, 0.02]
## 6.84e-04       | [0.00, 0.03]
## 0.06           | [0.01, 0.14]
# contrasts ta*fit effects between sessions (p-vals adjusted)
em4atr<-emtrends(m4a, pairwise ~ contingency, var = "ta")
em4atr <- summary(em4atr$contrasts)
em4atr
##  contrast          estimate    SE    df t.ratio p.value
##  (60/40) - (75/25)   0.0166 0.338 111.2   0.049  0.9987
##  (60/40) - (90/10)  -0.9694 0.368  77.5  -2.636  0.0271
##  (75/25) - (90/10)  -0.9860 0.334 110.5  -2.956  0.0106
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
effectsize::t_to_eta2(t=em4atr$t.ratio, df_error = em4atr$df, alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 2.17e-05       | [0.00, 0.01]
## 0.08           | [0.00, 0.22]
## 0.07           | [0.01, 0.18]

Plots

BIC scores (demeaned)

df = fit_data %>%
  group_by(ids,contingency,ta) %>%
  summarise_at(c("m1_BIC","m3_BIC"), mean, na.rm = TRUE)

# Demean 
df <- df %>% gather(model, BIC, m1_BIC:m3_BIC)
dem_df <- df %>%
  group_by(contingency) %>%
  summarise_at(c("BIC"), mean, na.rm = TRUE)
dem_df$BICmean <- dem_df$BIC
dem_df$BIC <- NaN
df<- merge(x = df, y = dem_df, by = "contingency", all.x = TRUE)
df$BIC_demean <- df$BIC.x - df$BICmean

df_summ = df %>%
  group_by(model, contingency) %>%
  summarise_each(funs(mean,sd), BIC_demean)

g <- ggplot(df_summ, aes(x=model, y=mean, fill=model)) + 
          geom_bar(inherit.aes = TRUE, stat = "summary", color="black", show.legend = FALSE, width=0.8, lwd=1) + 
          scale_fill_manual(values = pal2[c(15,16)]) +
          scale_x_discrete(labels=c("1-state", "n-state")) +
          labs(y= "BIC (demeaned)", x="")  +
          facet_grid(~contingency) +
          
          theme_bw() +
          raincloud_theme 
g
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

ggsave(here::here("output/figures/Fig_model_comparison.pdf"), width = 4, height = 4, dpi = 120)
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

Relative model fit and anxiety in 90/10

df <- fit_data[fit_data$contingency %in% c("90/10"),]
g <- ggplot(data = df, aes(x = m1m3_diff, y = ta_orig_scale, color=(m1m3_diff) )) +
geom_point( size = 5, alpha = 1, show.legend = FALSE) +
geom_smooth(method='lm', formula= y~x, alpha=0.3, show.legend = FALSE, color="black", fill="lightgray") +
#scale_color_gradient2(low=pal2[15], midpoint = 0, mid = "grey39", high=pal2[16]) +
  scale_color_gradient2(low=pal2[15], midpoint = 30, mid = pal2[16], high=pal2[16]) +
geom_vline(xintercept = 0, linetype = "longdash") +
theme_bw() +
stat_cor(method = "pearson",  alternative = "two.sided",  cor.coef.name = c("r"),  label.sep = ", ",  label.x.npc = "left",  label.y.npc = "top", digits = 3, r.digits = 3, p.digits = 3) + 
raincloud_theme +
labs(y= "Trait anxiety", x="Relative fit") 
g

ggsave(here::here("output/figures/Fig_model_fit_and_anxiety.pdf"), width = 7, height = 4.5, dpi = 120)


res<-cor.test(df$ta_orig_scale, df$m1m3_diff)
res
## 
##  Pearson's product-moment correlation
## 
## data:  df$ta_orig_scale and df$m1m3_diff
## t = 2.5135, df = 35, p-value = 0.01671
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.07673412 0.63464001
## sample estimates:
##       cor 
## 0.3910308
effectsize::t_to_eta2(t=res$statistic, df_error = res$parameter)
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.15           | [0.02, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

4b. Analysis of winning model parameters

N-state

data <- read.csv(here::here("data/models_no_states_and_switches_data.csv"))
 vars <- c("contingency", "ta", "study",  "ids")
 data <- assign_var_types(data, vars)
 df <- data %>% 
      melt(id.vars = c("ids", "ta", "study_str", "contingency"), measure.vars = c("tau_sh_reversalbsim", "tau_nosh_reversalbsim"), value.name = "tau", variable.name = "outcome" )
  df$outcome <- mapvalues(df$outcome,
                              from = c("tau_sh_reversalbsim", "tau_nosh_reversalbsim"),
                              to =c("sh", "nosh"))
  
m5a <- lmer(data=df, tau ~ contingency*outcome*ta + (1|ids) + (1|study_str) )
anova(m5a)
## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## contingency             0.3056  0.1528     2  92.534  0.5315    0.5895    
## outcome                10.6448 10.6448     1 228.265 37.0269 4.873e-09 ***
## ta                      0.0615  0.0615     1  82.771  0.2138    0.6450    
## contingency:outcome     1.1044  0.5522     2 228.265  1.9208    0.1488    
## contingency:ta          0.2530  0.1265     2 273.646  0.4400    0.6445    
## outcome:ta              0.0472  0.0472     1 228.265  0.1641    0.6858    
## contingency:outcome:ta  0.0249  0.0125     2 228.265  0.0434    0.9576    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 # effectsize
effectsize::effectsize(anova(m5a), type="eta", alternative="two.sided")
## # Effect Size for ANOVA (Type III)
## 
## Parameter              | Eta2 (partial) |       95% CI
## ------------------------------------------------------
## contingency            |           0.01 | [0.00, 0.07]
## outcome                |           0.14 | [0.07, 0.22]
## ta                     |       2.58e-03 | [0.00, 0.06]
## contingency:outcome    |           0.02 | [0.00, 0.06]
## contingency:ta         |       3.21e-03 | [0.00, 0.02]
## outcome:ta             |       7.18e-04 | [0.00, 0.02]
## contingency:outcome:ta |       3.80e-04 | [0.00, 0.00]
# statistical assumptions
## homogeneity of variance 
df$residuals <- residuals(m5a)
df$residuals_sq <- abs(residuals(m5a))^2
lev_m <- lm(residuals_sq ~ ids, data=df) #ANOVA of the squared residuals
anova(lev_m)
## Analysis of Variance Table
## 
## Response: residuals_sq
##            Df  Sum Sq  Mean Sq F value  Pr(>F)  
## ids        88  8.6164 0.097914  1.2494 0.09625 .
## Residuals 233 18.2597 0.078368                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## normality of residuals
ggplot(data=df, aes(x=residuals)) + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

car::qqPlot(df$residuals)

## [1] 70 66
 memm <- emmeans(m5a, specs = pairwise ~ outcome)
## NOTE: Results may be misleading due to involvement in interactions
 memm$emmeans
##  outcome emmean     SE   df lower.CL upper.CL
##  sh       1.125 0.0835 2.21    0.797     1.45
##  nosh     0.729 0.0835 2.21    0.401     1.06
## 
## Results are averaged over the levels of: contingency 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

1-state

data <- read.csv(here::here("data/models_no_states_and_switches_data.csv"))
vars <- c("contingency", "ta", "study",  "ids")
 data <- assign_var_types(data, vars)
 df <- data %>% 
      melt(id.vars = c("ids", "ta", "study_str", "contingency"), measure.vars = c("tau_sh_reversallbetav5", "tau_nosh_reversallbetav5"), value.name = "tau", variable.name = "outcome" )
  df$outcome <- mapvalues(df$outcome,
                              from = c("tau_sh_reversalbssm", "tau_nosh_reversalbssm"),
                              to =c("sh", "nosh"))
## The following `from` values were not present in `x`: tau_sh_reversalbssm, tau_nosh_reversalbssm
  m <- lmer(data=df, tau ~ contingency*outcome*ta + (1|ids) + (1|study_str) )
 anova(m)
## Type III Analysis of Variance Table with Satterthwaite's method
##                        Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## contingency            1.3266  0.6633     2 164.13  2.5281   0.08291 .  
## outcome                6.5164  6.5164     1 227.86 24.8369 1.234e-06 ***
## ta                     0.1749  0.1749     1  91.99  0.6665   0.41639    
## contingency:outcome    0.6815  0.3408     2 227.86  1.2988   0.27488    
## contingency:ta         0.3015  0.1508     2 265.93  0.5746   0.56363    
## outcome:ta             0.0854  0.0854     1 227.86  0.3256   0.56882    
## contingency:outcome:ta 0.0698  0.0349     2 227.86  0.1330   0.87552    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 # effectsize
effectsize::effectsize(anova(m), type="eta", alternative="two.sided")
## # Effect Size for ANOVA (Type III)
## 
## Parameter              | Eta2 (partial) |       95% CI
## ------------------------------------------------------
## contingency            |           0.03 | [0.00, 0.09]
## outcome                |           0.10 | [0.04, 0.18]
## ta                     |       7.19e-03 | [0.00, 0.08]
## contingency:outcome    |           0.01 | [0.00, 0.05]
## contingency:ta         |       4.30e-03 | [0.00, 0.03]
## outcome:ta             |       1.43e-03 | [0.00, 0.03]
## contingency:outcome:ta |       1.17e-03 | [0.00, 0.01]
 memm <- emmeans(m, specs = pairwise ~ outcome)
## NOTE: Results may be misleading due to involvement in interactions
 memm$emmeans
##  outcome                  emmean     SE   df lower.CL upper.CL
##  tau_sh_reversallbetav5    1.166 0.0915 2.44    0.832     1.50
##  tau_nosh_reversallbetav5  0.856 0.0915 2.44    0.522     1.19
## 
## Results are averaged over the levels of: contingency 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
 memm <- emmeans(m, specs = pairwise ~ outcome*contingency)
## NOTE: Results may be misleading due to involvement in interactions
 memm$emmeans
##  outcome                  contingency emmean     SE   df lower.CL upper.CL
##  tau_sh_reversallbetav5   60/40        1.037 0.1255 7.60    0.744     1.33
##  tau_nosh_reversallbetav5 60/40        0.857 0.1255 7.60    0.565     1.15
##  tau_sh_reversallbetav5   75/25        1.172 0.0866 3.01    0.897     1.45
##  tau_nosh_reversallbetav5 75/25        0.761 0.0866 3.01    0.486     1.04
##  tau_sh_reversallbetav5   90/10        1.289 0.1244 7.31    0.997     1.58
##  tau_nosh_reversallbetav5 90/10        0.950 0.1244 7.31    0.658     1.24
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

5. Analysis of model-free learning rates from oddball and meaningful events

Stats

Oddball/meaninful by anx and outcome

# Load calcualted model-free learning rates
mflr_data <- read.csv(here::here("data/oddball_data_n5.csv"))
mflr_data = assign_var_types(mflr_data, c("study", "lrtype", "ta","contingency", "ids" ))



df <- melt(mflr_data, measure.vars=c("mflr_sh","mflr_nosh"),
    value.name = "mflr", variable.name="outcome_str")
df = assign_var_types(df, c("outcome_str" ))


## LOG learningr ates for analysis

df$log_mflr <- log(df$mflr)
df_stats<-df[!is.na(df$log_mflr) & !is.infinite(df$log_mflr),]

m6a<-lmer(log_mflr ~ ta*contingency*lrtype*outcome_str + (1|ids) + (1|study_str), df_stats)
## boundary (singular) fit: see help('isSingular')
anova(m6a)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                    Sum Sq Mean Sq NumDF  DenDF F value
## ta                                 0.0001  0.0001     1  97.48  0.0004
## contingency                        1.7758  0.8879     2 587.42  2.2837
## lrtype                             7.0159  7.0159     1 533.22 18.0453
## outcome_str                       16.0828 16.0828     1 533.30 41.3657
## ta:contingency                     0.8964  0.4482     2 589.62  1.1527
## ta:lrtype                          3.2520  3.2520     1 533.52  8.3643
## contingency:lrtype                 1.9219  0.9609     2 533.20  2.4716
## ta:outcome_str                     0.9197  0.9197     1 533.67  2.3654
## contingency:outcome_str            1.7216  0.8608     2 533.21  2.2140
## lrtype:outcome_str                 0.0794  0.0794     1 533.30  0.2043
## ta:contingency:lrtype              1.6488  0.8244     2 533.30  2.1204
## ta:contingency:outcome_str         1.5256  0.7628     2 533.25  1.9620
## ta:lrtype:outcome_str              0.3092  0.3092     1 533.67  0.7954
## contingency:lrtype:outcome_str     0.0866  0.0433     2 533.21  0.1114
## ta:contingency:lrtype:outcome_str  0.1033  0.0517     2 533.25  0.1328
##                                      Pr(>F)    
## ta                                 0.984937    
## contingency                        0.102810    
## lrtype                            2.546e-05 ***
## outcome_str                       2.810e-10 ***
## ta:contingency                     0.316483    
## ta:lrtype                          0.003983 ** 
## contingency:lrtype                 0.085416 .  
## ta:outcome_str                     0.124645    
## contingency:outcome_str            0.110262    
## lrtype:outcome_str                 0.651486    
## ta:contingency:lrtype              0.120988    
## ta:contingency:outcome_str         0.141589    
## ta:lrtype:outcome_str              0.372888    
## contingency:lrtype:outcome_str     0.894576    
## ta:contingency:lrtype:outcome_str  0.875626    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effectsize
effectsize::effectsize(anova(m6a), type="eta", alternative="two.sided")
## # Effect Size for ANOVA (Type III)
## 
## Parameter                         | Eta2 (partial) |       95% CI
## -----------------------------------------------------------------
## ta                                |       3.68e-06 | [0.00, 0.00]
## contingency                       |       7.72e-03 | [0.00, 0.03]
## lrtype                            |           0.03 | [0.01, 0.07]
## outcome_str                       |           0.07 | [0.04, 0.12]
## ta:contingency                    |       3.89e-03 | [0.00, 0.02]
## ta:lrtype                         |           0.02 | [0.00, 0.04]
## contingency:lrtype                |       9.19e-03 | [0.00, 0.03]
## ta:outcome_str                    |       4.41e-03 | [0.00, 0.02]
## contingency:outcome_str           |       8.24e-03 | [0.00, 0.03]
## lrtype:outcome_str                |       3.83e-04 | [0.00, 0.01]
## ta:contingency:lrtype             |       7.89e-03 | [0.00, 0.03]
## ta:contingency:outcome_str        |       7.30e-03 | [0.00, 0.03]
## ta:lrtype:outcome_str             |       1.49e-03 | [0.00, 0.01]
## contingency:lrtype:outcome_str    |       4.18e-04 | [0.00, 0.01]
## ta:contingency:lrtype:outcome_str |       4.98e-04 | [0.00, 0.01]
# statistical assumptions
## homogeneity of variance 
df_stats$residuals <- residuals(m6a)
df_stats$residuals_sq <- abs(residuals(m6a))^2
lev_m <- lm(residuals_sq ~ ids, data=df_stats) #ANOVA of the squared residuals
anova(lev_m)
## Analysis of Variance Table
## 
## Response: residuals_sq
##            Df Sum Sq Mean Sq F value   Pr(>F)   
## ids        88  45.37 0.51556  1.4306 0.009747 **
## Residuals 544 196.05 0.36038                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
levRes = car::leveneTest(residuals_sq ~ ids, data=df_stats)
levRes
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group  88  1.4094 0.01266 *
##       544                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## normality of residuals
ggplot(data=df_stats, aes(x=residuals)) + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

car::qqPlot(df_stats$residuals)

## [1] 236  80
# post hocs
em5a1 = emmeans(m6a, specs = pairwise ~ contingency:outcome_str)
## NOTE: Results may be misleading due to involvement in interactions
em5a1$emmeans
##  contingency outcome_str emmean     SE   df lower.CL upper.CL
##  60/40       mflr_sh      -1.49 0.1035 7.45    -1.73    -1.25
##  75/25       mflr_sh      -1.40 0.0603 3.92    -1.57    -1.23
##  90/10       mflr_sh      -1.52 0.1033 7.38    -1.76    -1.28
##  60/40       mflr_nosh    -1.67 0.1035 7.45    -1.91    -1.43
##  75/25       mflr_nosh    -1.84 0.0606 3.97    -2.01    -1.67
##  90/10       mflr_nosh    -1.94 0.1028 7.23    -2.18    -1.70
## 
## Results are averaged over the levels of: lrtype 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
em5a1$contrasts
##  contrast                              estimate     SE  df t.ratio p.value
##  (60/40 mflr_sh) - (75/25 mflr_sh)      -0.0876 0.1037 130  -0.845  0.9585
##  (60/40 mflr_sh) - (90/10 mflr_sh)       0.0291 0.1043 527   0.279  0.9998
##  (60/40 mflr_sh) - (60/40 mflr_nosh)     0.1829 0.1042 525   1.755  0.4963
##  (60/40 mflr_sh) - (75/25 mflr_nosh)     0.3473 0.1036 135   3.353  0.0130
##  (60/40 mflr_sh) - (90/10 mflr_nosh)     0.4530 0.1039 527   4.360  0.0002
##  (75/25 mflr_sh) - (90/10 mflr_sh)       0.1166 0.1034 130   1.128  0.8691
##  (75/25 mflr_sh) - (60/40 mflr_nosh)     0.2704 0.1037 130   2.609  0.1024
##  (75/25 mflr_sh) - (75/25 mflr_nosh)     0.4348 0.0674 526   6.451  <.0001
##  (75/25 mflr_sh) - (90/10 mflr_nosh)     0.5406 0.1030 128   5.250  <.0001
##  (90/10 mflr_sh) - (60/40 mflr_nosh)     0.1538 0.1043 527   1.474  0.6808
##  (90/10 mflr_sh) - (75/25 mflr_nosh)     0.3182 0.1033 134   3.080  0.0296
##  (90/10 mflr_sh) - (90/10 mflr_nosh)     0.4240 0.1037 526   4.087  0.0007
##  (60/40 mflr_nosh) - (75/25 mflr_nosh)   0.1644 0.1036 135   1.588  0.6081
##  (60/40 mflr_nosh) - (90/10 mflr_nosh)   0.2702 0.1039 527   2.600  0.0990
##  (75/25 mflr_nosh) - (90/10 mflr_nosh)   0.1058 0.1028 132   1.028  0.9076
## 
## Results are averaged over the levels of: lrtype 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 6 estimates
em5a2 = emmeans(m6a, specs = pairwise ~ outcome_str)
## NOTE: Results may be misleading due to involvement in interactions
em5a2$emmeans
##  outcome_str emmean     SE   df lower.CL upper.CL
##  mflr_sh      -1.47 0.0689 1.78    -1.80    -1.14
##  mflr_nosh    -1.82 0.0690 1.78    -2.15    -1.48
## 
## Results are averaged over the levels of: contingency, lrtype 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
em5a2$contrasts
##  contrast            estimate     SE  df t.ratio p.value
##  mflr_sh - mflr_nosh    0.347 0.0539 526   6.439  <.0001
## 
## Results are averaged over the levels of: contingency, lrtype 
## Degrees-of-freedom method: kenward-roger
# relationship with anxiety
emtrm6<-emtrends(m6a, pairwise ~ lrtype, var = "ta")
## NOTE: Results may be misleading due to involvement in interactions
emtrm6 <- summary(emtrm6$contrasts)
emtrm6
##  contrast     estimate      SE  df t.ratio p.value
##  common - odd   0.0141 0.00487 526   2.892  0.0040
## 
## Results are averaged over the levels of: contingency, outcome_str 
## Degrees-of-freedom method: kenward-roger
effectsize::t_to_eta2(t=emtrm6$t.ratio, df_error = emtrm6$df, alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.02           | [0.00, 0.04]

Oddball/meaningful by best model per participant (LOG)

data <- read.csv(here::here("data/oddball_data_n5.csv"))
data$mflr_overall_log <- log(data$mflr_overall) 
data$lrtype = to_factor(data$lrtype)
data$ids = to_factor(data$ids)
data$ta <- demean(data$ta)
data["tabin"]<-dicho(data$ta)
data$tabin <- mapvalues(data$tabin,
                           from = c(0,1),
                           to = c("lowAnx", "hiAnx"))

odd<-data[data$lrtype %in% "odd",]
common<-data[data$lrtype %in% "common",]
data <- merge(odd, common, by="specID", all.x = TRUE)

data$mflr_diff <- (data$mflr_overall.y) - (data$mflr_overall.x)
data$mflr_diff_log <- (data$mflr_overall_log.y) - (data$mflr_overall_log.x)
data<-data[!is.na(data$mflr_diff_log) & !is.infinite(data$mflr_diff_log),]

data <- data[, -grep(".y", colnames(data))]
colnames(data)<-gsub(".x","",colnames(data))
data <- data %>%
  group_by(ids) %>%
  summarise_at(c("mflr_diff_log", "mflr_diff"), mean, na.rm = TRUE)


fitd<-read.csv(here::here("data/model_fit_final_full.csv"))
fitd<-assign_var_types(fitd, c("contingency", "ta", "study", "ids"))   
fitd <- fitd %>%
  group_by(ids, ta, study_str) %>%
  summarise_at(c("m1_BIC", "m3_BIC", "m1m3_diff"), mean, na.rm = TRUE)

fitd["best_model"] <- apply(fitd[,c("m1_BIC", "m3_BIC")], 1, which.min)
fitd$best_model <- mapvalues(fitd$best_model,
                            from = c(1,2),
                            to =c("1-state", "n-state"))
fitd <- fitd[,c("ids", "best_model", "study_str", "m1m3_diff")]
data <-merge(data, fitd, by="ids")

## analysis
data$best_model <- as.factor(data$best_model)
tt = t.test(mflr_diff_log~best_model, alternative = "two.sided", var.equal = FALSE, data=data)
tt
## 
##  Welch Two Sample t-test
## 
## data:  mflr_diff_log by best_model
## t = -2.107, df = 68.788, p-value = 0.03877
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.321982745 -0.008787365
## sample estimates:
## mean in group 1-state mean in group n-state 
##             0.1295605             0.2949456
effectsize::t_to_eta2(t=tt$statistic, df_error = tt$parameter, alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.06           | [0.00, 0.20]
 data %>%
  group_by(best_model) %>%
  summarise_at(c("mflr_diff_log", "mflr_diff"), mean, na.rm = TRUE)
## # A tibble: 2 × 3
##   best_model mflr_diff_log mflr_diff
##   <fct>              <dbl>     <dbl>
## 1 1-state            0.130    0.0212
## 2 n-state            0.295    0.0586
cor.test(data$mflr_diff_log, data$m1m3_diff, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  data$mflr_diff_log and data$m1m3_diff
## S = 90596, p-value = 0.03121
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.228839
confintr::ci_cor(data$mflr_diff_log, data$m1m3_diff, method = "spearman",type="bootstrap")
## 
##  Two-sided 95% bootstrap confidence interval for the true Spearman
##  correlation coefficient based on 9999 bootstrap replications and the
##  bca method
## 
## Sample estimate: 0.228839 
## Confidence interval:
##        2.5%       97.5% 
## 0.005350783 0.428545212
odd_model_data = data
write.csv(odd_model_data, here::here("data",  "oddball_by_model.csv"))

Oddball/meaningful by anxiety (LOG)

df2 <- df %>%
  group_by(tabin, lrtype, ids, ta_orig_scale) %>%
  summarise_at(c("log_mflr"), mean, na.rm = TRUE)

dfnew <- df2 %>% tidyr::pivot_wider(
                  id_cols=c("ids", , "tabin", "ta_orig_scale"),
                  names_from=lrtype,
                  values_from=log_mflr 
                  ) 
dfnew$mflr_diff <- dfnew$common - dfnew$odd
dfnew<-dfnew[!is.na(dfnew$mflr_diff) & !is.infinite(dfnew$mflr_diff),]

# somewhat artificially also do t-test
tt = t.test(mflr_diff~tabin, alternative = "two.sided", var.equal = FALSE, data=dfnew)
tt
## 
##  Welch Two Sample t-test
## 
## data:  mflr_diff by tabin
## t = -2.2445, df = 77.049, p-value = 0.02767
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.34770972 -0.02078966
## sample estimates:
## mean in group lowAnx  mean in group hiAnx 
##            0.1691060            0.3533557
effectsize::t_to_eta2(t=tt$statistic, df_error = tt$parameter, alternative = "two.sided")
## Eta2 (partial) |       95% CI
## -----------------------------
## 0.06           | [0.00, 0.19]
dfnew %>%
  group_by(tabin) %>%
  summarise_at(c("mflr_diff"), mean, na.rm = TRUE)
## # A tibble: 2 × 2
##   tabin  mflr_diff
##   <fct>      <dbl>
## 1 lowAnx     0.169
## 2 hiAnx      0.353
cor.test(dfnew$mflr_diff, dfnew$ta_orig_scale, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  dfnew$mflr_diff and dfnew$ta_orig_scale
## S = 81629, p-value = 0.03324
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2298771
confintr::ci_cor(dfnew$mflr_diff, dfnew$ta_orig_scale, method = "spearman",type="bootstrap")
## 
##  Two-sided 95% bootstrap confidence interval for the true Spearman
##  correlation coefficient based on 9999 bootstrap replications and the
##  bca method
## 
## Sample estimate: 0.2298771 
## Confidence interval:
##       2.5%      97.5% 
## 0.03696545 0.41365519

Plots

Oddball/Meaninful by anxiety

df2 <- df %>%
  group_by(tabin, lrtype, ids, ta_orig_scale) %>%
  summarise_at(c("mflr"), mean, na.rm = TRUE)

dfnew <- df2 %>% tidyr::pivot_wider(
                  id_cols=c("ids", , "tabin", "ta_orig_scale"),
                  names_from=lrtype,
                  values_from=mflr 
                  ) 
dfnew$mflr_diff <- dfnew$common - dfnew$odd

g <- ggplot(data = dfnew, aes(y = mflr_diff, x = tabin, fill = tabin)) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2, show.legend = TRUE) +
geom_point(aes(color= tabin), position =   position_jitterdodge(  jitter.width = .15,  dodge.width = 0.25), size = 2, alpha = 0.5, show.legend=FALSE) + 
geom_boxplot(width = .25, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=5, shape=23, aes(group=tabin), color="black", 
             position = position_dodge( 0.25), stroke=1, show.legend = FALSE) + 
expand_limits(x = 3) +
#scale_color_manual(values = pal3[c(2,4)] , name = "", labels = c("Low Anx", "High Anx")) +
#scale_fill_manual(values = pal3[c(2,4)], name = "", labels = c("Low Anx", "High Anx")) + 
scale_color_manual(values = c("olivedrab2", "darkgreen") , name = "", labels = c("Low Anx", "High Anx")) +
scale_fill_manual(values = c("olivedrab2", "darkgreen"), name = "", labels = c("Low Anx", "High Anx")) + 
theme_bw() +
raincloud_theme +
  theme(strip.text.x = element_text(size = 14,  face = "bold"),
      strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
theme(legend.position = "top") + 
labs(y= "Meaningful-oddball learning difference", x="Trait anxiety") +
  scale_x_discrete(labels=c("Low", "High")) +
guides(fill=guide_legend(nrow=4,byrow=TRUE)) + 
  ylim(-0.1, 0.35)
g

ggsave(here::here("output/figures/Fig_oddball_anxiety.pdf"), width = 4, height = 5, dpi = 120)

Also plot continously

df2 <- df %>%
  group_by(tabin, lrtype, ids, ta_orig_scale) %>%
  summarise_at(c("log_mflr"), mean, na.rm = TRUE)

dfnew <- df2 %>% tidyr::pivot_wider(
                  id_cols=c("ids", , "tabin", "ta_orig_scale"),
                  names_from=lrtype,
                  values_from=log_mflr 
                  ) 
dfnew$mflr_diff <- dfnew$common - dfnew$odd
dfnew<-dfnew[!is.na(dfnew$mflr_diff) & !is.infinite(dfnew$mflr_diff),]

g <- ggplot(data = dfnew, aes(x = ta_orig_scale, y = mflr_diff, color=ta_orig_scale ), show.legend = FALSE) +
geom_point(size=3 ) +
geom_smooth(method='lm', formula= y~x, alpha=0.3, show.legend = FALSE, color="black") +
scale_color_gradient2(low="olivedrab2", high="darkgreen", midpoint = 40, mid = "olivedrab") +
  #scale_color_gradient2(low=pal2[15], midpoint = 0, mid = "grey39", high=pal2[16]) +
theme_bw() +
raincloud_theme +
  stat_cor(method = "spearman",  alternative = "two.sided",  cor.coef.name = c("r"),  label.sep = ", ",  label.x.npc = "left",  label.y.npc = "top", digits = 3, r.digits = 3, p.digits = 3) + 
  
  labs(x= "Trait anxiety", y="Meaningful-oddball learning difference")   
g

ggsave(here::here("output/figures/Fig_oddball_anxiety_continuous.pdf"), width = 6.5, height = 5, dpi = 120)
df2 <- df %>%
  group_by(contingency, tabin, lrtype, ids) %>%
  summarise_at(c("log_mflr"), mean, na.rm = TRUE)



g <- ggplot(data = df2, aes(y = log_mflr, x = tabin, fill = interaction(tabin, lrtype))) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2, show.legend = TRUE) +
geom_point(aes(color= interaction( tabin, lrtype)), position =   position_jitterdodge(  jitter.width = .15,  dodge.width = 0.25), size = 2, alpha = 0.5, show.legend=FALSE) + 
geom_boxplot(width = .2, outlier.shape = NA, alpha = 0.7, lwd=1.2, show.legend = FALSE) +
stat_summary(geom="point", fun="mean", size=5, shape=23, aes(group= interaction(tabin, lrtype)), color="black", 
             position = position_dodge( 0.25), stroke=1, show.legend = FALSE) + 
expand_limits(x = 3) +
scale_color_manual(values = pal[c(2,4,1,3,2,4)] ) +
scale_fill_manual(values = pal[c(2,4,1,3,2,4)]) + 
theme_bw() +
raincloud_theme +
theme(legend.position = "right") + 
labs(y= "Learning rate", x="") +
  scale_x_discrete(labels=c("common" = "Meaningful", "odd" =  "Oddball")) + 
guides(fill=guide_legend(nrow=4,byrow=TRUE)) 
g

Oddball/Meaninful by best fitted model

## Plot the difference as well
g <- ggplot(data = odd_model_data, aes(y = mflr_diff, x = best_model, fill =  best_model)) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8, lwd=1.2) +
geom_point(aes(color= best_model), position =   position_jitterdodge(  jitter.width = .45,  dodge.width = 1), size = 2, alpha = 0.5, show.legend=FALSE) +  
geom_boxplot(inherit.aes = TRUE, width = 0.25, lwd=1.2, outlier.shape = NA, alpha = 0.7, show.legend = FALSE) + 
stat_summary(geom="point", fun="mean", size=5, shape=23, aes(group=best_model), color="black", 
             position = position_dodge( 1), stroke=1, show.legend = FALSE) + 
expand_limits(x = 3) +
scale_color_manual(values = c(pal2[15], pal2[16], pal[4], pal[3])) +
scale_fill_manual(values = c(pal2[15], pal2[16], pal[4], pal[3]), name = "", labels = c("1-state", "n-state")) +
theme_bw() +
raincloud_theme +
theme(strip.text.x = element_text(size = 14,  face = "bold"),
      strip.background = element_rect(color="black", fill="#ebebeb", size=0, linetype="solid")) +
theme(legend.position = "top") +  #c(0.8, 0.2)

labs(y= "Meaningful-oddball learning difference", x="Best model")  +

  guides(fill=guide_legend(nrow=1,byrow=TRUE)) +
ylim(-0.1, 0.35)
g

ggsave(here::here("output/figures/Fig_oddball_model.pdf"), width = 4, height = 5, dpi = 120)
g <- ggplot(data = odd_model_data, aes(x = m1m3_diff, y = mflr_diff_log, color=m1m3_diff ), show.legend = FALSE) +
geom_point(size=3 ) +
geom_smooth(method='lm', formula= y~x, alpha=0.3, show.legend = FALSE, color="black") +
#scale_color_gradient2(low=pal[3], high=pal[4], midpoint = -5, mid = pal[3]) +
  scale_color_gradient2(low=pal2[15], midpoint = 30, mid = pal2[16], high=pal2[16]) +
  #scale_color_gradient2(low=pal2[15], midpoint = 0, mid = "grey39", high=pal2[16]) +
   geom_vline(xintercept = 0, linetype = "longdash") +
theme_bw() +
raincloud_theme +
  stat_cor(method = "spearman",  alternative = "two.sided",  cor.coef.name = c("r"),  label.sep = ", ",  label.x.npc = "left",  label.y.npc = "top", digits = 3, r.digits = 3, p.digits = 3) + 
  labs(x= "Model difference", y="Meaningful-oddball learning difference")   
g

ggsave(here::here("output/figures/Fig_oddball_model_continuous.pdf"), width = 6.5, height = 5, dpi = 120)